The statistics of identifying differentially expressed genes in Expresso and TM4: a comparison
© Sioson et al; licensee BioMed Central Ltd. 2006
Received: 16 August 2005
Accepted: 20 April 2006
Published: 20 April 2006
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.
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.
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.
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 . One such method is the two-stage, log-linear ANOVA mixed models technique of Wolfinger, et al., . 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  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  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.,  compare two log-linear 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  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.,  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 ), naive methods, such as total intensity normalization, are still in use (e.g., Held, et al., ). Xie, et al.,  did a comparative study of normalization methods and test statistics to analyze the results of a DNA-protein binding microarray experiment. Using performance and bias correction criteria, Bolstad, et al.,  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.,  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.,  describe the features and capabilities of TM4, while Quackenbush  describes the normalization and transformation methods implemented in it. Williams, et al., , Zhu, et al., , and Khaitovich, et al.,  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–20]. Expresso has adapted the two-stage ANOVA mixed models technique of Wolfinger, et al.,  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 . 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
Normalization and low intensity filtering in TM4
Quackenbush  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 IIVTL 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.
Number of differentially expressed genes. Numbers of identified differentially expressed genes in the GP and GOT models after preprocessing by 0 or more MIDAS pipeline steps — IIV, IIVT, IIVTL, IIVTLS, or IIVTSLF. 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).
GP Model — Experiment 1
GOT Model — Experiment 2
Retention counts and percentages. Retention counts (RC) and retention percentages (RP) for the differentially expressed gene sets of Table 1. RC is the number of genes in a set before the MIDAS step that remain in the set after that step. RC is the number of genes in the set before the MIDAS step that remain in the set after that step. RP is the percentage of remaining genes with respect to the number of genes in the set after the MIDAS step. RC and RP are reported for the intersections IIV∩ IIVT, IIVT ∩ IIVTL, IIVTL ∩ IIVTLS, and IIVTLS ∩ IIVTLSF, as well as intersection IIV∩ ∩IIVTLSF, which corresponds to the effect of MIDAS steps from the start of the pipeline to the end. The ALL entries correspond to the H1 hypotheses of the GOT model (see Materials and Methods).
GP Model — Experiment 1
GOT Model — Experiment 2
In Experiment 1 results, the number of genes commonly assessed by Expresso as significantly expressed when using IIV and IIVT is high. For example, there is 95.45% retention of WT genes (545 total) assessed as up-expressed when using IIVT 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 IIVTS data in Expresso compared to when using IIVT. While we observe increase in the retention percentage in IIVTLS (from IIVTL) and IIVTLSF (from IIVTLS), there's low retention percentage in the results using IIVTLSF from IIV data. This can be traced in the low retention percentage of IIVTL from IIVT. 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 ∩ IIVTLSF 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.
Comparison of MIV and IIV. We compare MIV and IIV input data with respect to the sets of genes ultimately identified as differentially expressed. The genotype and set labelings are the same as those in Table 1. The GP model was used to analyze the unnormalized MIV and IIV data from Experiment 1. The GOT model was used to analyze the unnormalized MIV and IIV data from Experiment 2. Up-expressed and down-expressed counts are reported for both experiments and all genotypes. The count of genes in the intersections is found in column Common. For convenience, percentages of the intersection with respect to the MIV and IIV sets are tabulated in the last two columns. The ALL entries correspond to the H1 hypotheses of the GOT model (see Materials and Methods).
% of common in MIV
% of common in IIV
GP Model — Experiment 1
GOT Model — Experiment 2
Comparison of statistical methods
Comparison of Expresso analysis and MEV t test. We compare the number of genes identified as differentially expressed by Expresso analysis and the MEV t-test. Counts for both up-expressed and down-expressed genes, as well as all four genotypes of Experiment 2, are reported. The first three analyses — the GP model, the GOT model, and the MEV t-test — take the IIVTLSF data as input. For point of comparison, the last analysis uses the GP model on MIV input.
Comparison of qRT-PCR with Expresso and TM4. We compare the qRT-PCR results with identified up-expressed and down-expressed genes in Experiment 2 using Expresso and TM4. Results of qRT-PCR are available for each of the four ecotypes in the numbers n in the second column. Numbers of agreement or non-agreement are shown in the S, D, and F columns. The S (same) column tabulates the number of genes for which the sign of the log(fold change) for statistically significant differential expression matches the direction of change in the corresponding qRT-PCR result. The D (differing) column tabulates the number of genes for which the sign of the log(fold change) for statistically significant differential expression is in the opposite direction of the change in the corresponding qRT-PCR result. The F (filtered) column tabulates the number of genes for which there is a qRT-PCR result, but for which either the gene was filtered by MIDAS low intensity filtering or the analysis method did not assess the change in expression as statistically significant. The MEV t-test (first grouping) results are for the typical TM4 process, which involves IIV input data followed by the four MIDAS steps, which we denote IIVTLSF The GP model (second grouping) gives the same numbers and uses the same IIVTLSFinput data, the GP model (third grouping) results use MIV input data and has no filtered genes.
Correlation of qRT-PCR with Expresso and TM4. We calculate the correlation of qRT-PCR results from the Expresso GP model and the MEV t-test. For each comparison, the analytical result for each gene for which there is a qRT-PCR result for a particular genotype is assembled into a result vector. We used SAS to compute a Pearson correlation of each result vector with the corresponding vector of qRT-PCR results. The computed correlations are as reported above.
qRT-PCR versus GP model with MIV input
qRT-PCR versus GP model with IIVTLSF input
qRT-PCR versus MEV t-test with IIVTLSF input
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 intensity-dependent 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 (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.
Microarray data sets
We used data sets from two experiments that used the Arabidopsis Oligonucleotide Microarrays , 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.
Experiment 1 compares the responses of Arabidopsis thaliana wild type, ecotype Columbia (henceforth, WT), and of an antisense plant for phospholipase D α (antiPLD) in Columbia background to drought stress . Plants were harvested at a single time point, and two biological replicate hybridizations were done for each of WT and antiPLD. Scan Array Express (PerkinElmer Life and Analytical Sciences, Inc., Boston, MA USA) was used to quantitate the four microarrays. By default, ScanArray Express performs a global lowess normalization of median intensities per microarray. ScanArray Express output its results to four files in GenePix (GPR) format, which constitute the Experiment 1 data set.
Li, et al.,  compare the responses to elevated CO2 of a wild Arabidopsis thaliana relative (Thellungiella halophila, ecotype Shandong; Th) and of three Arabidopsis thaliana ecotypes: Wassilewskija (WS), Columbia (Col-0), and Cape Verde Islands (Cvi-0). Three biological replicate hybridizations were done for each genotype. GenePix (Axon Instruments, Union City, CA USA) was used to quantitate the twelve microarrays. GenePix also performs by default a global lowess normalization of median intensities per microarray. The output of GenePix is twelve GPR files, which constitute the Experiment 2 data set.
Real-time quantitative RT-PCR
For verification of microarray results in Experiment 2, Li, et al.,  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., . 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.,  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.
Expresso GP model
The normalization model is
y spdab = μ + P p + D d + A a + (P × A) pa + B ba + r spdab .
Each y spdab value is the log2-transformed intensity of spot s within the dye d image in block b of array a. (A spot may represent a gene, a control, or a blank.) The global mean of the y spdab values, over all microarrays, is μ. The fixed effects in the model are the plant sample effect P p , where p indexes the various distinct plant samples from which mRNA was obtained, and the dye effect D d , where d has two values for the two dyes. The random effects in the model are the array effect A a , where a indexes the microarrays, the interaction effect (P × A) pa of plant sample p with microarray a, and the block effect B ba , where b identifies the block within microarray a. The model residual is r spdab . This differs from the normalization model in Wolfinger, et al., , in that it incorporates dye and block effects. It is a refinement of the Expresso normalization model in Watkinson, et al., , 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
r gspdab = G g + (G × P) gp + (G × D) gd + (G × A) ga + λ gspdab .
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., , in that it incorporates interactions between gene and dye and between gene and array. It refines the Expresso gene model in Watkinson, et al.,  to include the interaction between gene and array.
The estimate of the expression level of each gene in each treatment comparison is done by computing the pair-wise least square mean differences of gene-treatment effects. The Tukey-Kramer multiple comparison of gene-(plant sample) effects on each gene is made to estimate the p values associated with each calculated expression level. If there are ρ plant samples, then there are possible pairwise comparisons. If we index the plant samples from 1 to ρ, then the null hypothesis for gene g and comparison i, j, where 1 ≤ i <j ≤ ρ, is
Ho: (G × P) gi = (G × P) gi .
The difference (G × P) gi - (G × P) gj is the estimate of the log2(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.
The above GP model was used to analyze both Experiment 1 and Experiment 2. In both experiments, there are 48 blocks per array. In Experiment 1, there are four arrays and four plant samples, namely, WT-control, WT-stressed, antiPLD-control, and antiPLD-stressed. In Experiment 2, there are 12 arrays and eight plant samples, namely, Col-0-test, Col-0-control, Cvi-0-test, Cvi-0-control, WS-test, WS-control, Th-test, and Th-control.
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 gene-genotype-treatment model, an alternative set of log-linear ANOVA mixed models, for the elevated CO2 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
y sotdab = μ + O o + TT + D d + (O × T) ot + B ba + r sotdab .
Each y sotdab value is the log2-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. The fixed effects in the model are the genotype effect O o , where o indexes the genotype (organism), the treatment effect T t , where t is the treatment, and the dye effect D d , where d has two values for the two dyes. The random effects in the model are the interaction effect (O × T) ot of genotype o with microarray a, and the block effect B ba , where b identifies the block within microarray a. The model residual is r sotdab . This differs from the normalization model in Wolfinger, et al., , in that it incorporates genotype (organism), dye, and block effects. It is a refinement of the Expresso normalization model in Watkinson, et al., , 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
r gsotdab = G g + (G × O) go + (G × T) gt + (G × O × T) got + (G × D) gd + λ gsotdab
Here, G g are as for the GP model. The interactions (G × O) go of gene g with genotype o, (G × T) gt of gene g with treatment t, (G × O × T) got of gene g with genotype o and treatment t, and (G × D) gd of gene g with dye d are fixed effects. The λ gsotdab values are stochastic errors. This differs from the gene model in Wolfinger, et al., , 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.,  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
H1: (G × T) gi = (G × T) gj .
The difference (G × T) gi - (G × T) gj is the estimate of the log2(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
H2: (G × O × T) goi = (G × O × T) goj .
The difference (G × O × T) goi - (G × O × T) goj is the estimate of the log2 (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 . 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.
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  discusses this normalization in the context of sev-eral variations that are possible to address differing channel intensities.
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 = log10(R i G i ) and y i = log2(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  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 .
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 per-block variances of log(R i /G i ) values are the same [25, 28]. Quackenbush  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  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.
The work has been supported by the National Science Foundation (DBI-0223905 and BIO/IBN 0219322) and Virginia Tech and UIUC institutional funds.
- Churchill GA: Using ANOVA to Analyze Microarray Data. Bio Techniques 2004, 37(2):173–177.Google Scholar
- Wolfinger R, Gibson G, Wolfinger E, Bennett L, Hamadeh H, Bushel P, Afshari C, Paules R: Assessing Gene Significance from cDNA Microarray Expression Data via Mixed Models. Journal of Computational Biology 2001, 8: 625–637. 10.1089/106652701753307520View ArticlePubMedGoogle Scholar
- Kerr MK: Linear Models for Microarray Data Analysis: Hidden Similarities and Differences. Journal of Computational Biology 2003, 10(6):891–901. 10.1089/106652703322756131View ArticlePubMedGoogle Scholar
- Pan W: A Comparative Review of Statistical Methods for Discovering Differentially Expressed Genes in Replicated Microarray Experiments. Bioinformatics 2002, 18(4):546–554. 10.1093/bioinformatics/18.4.546View ArticlePubMedGoogle Scholar
- Chu TM, Weir B, Wolfinger RD: Comparison of Li-Wong and Loglinear Mixed Models for the Statistical Analysis of Oligonucleotide Arrays. Bioinformatics 2004, 20(4):500–506. 10.1093/bioinformatics/btg435View ArticlePubMedGoogle Scholar
- Cui X, Churchill GA: Statistical Tests for Differential Expression in cDNA Microarray Experiments. Genome Biology 2003., 4(210):
- Rosa GJ, Steibel JP, Tempelman RJ: Reassessing Design and Analysis of Two-colour Microarray Experiments using Mixed Effects Models. Comparative and Functional Genomics 2005, 6: 123–131. 10.1002/cfg.464PubMed CentralView ArticlePubMedGoogle Scholar
- Futschik M, Crompton T: Model Selection and Efficiency Testing for Normalization of cDNA Microarray Data. Genome Biology 2004., 5(R60):
- Held M, Gase K, Baldwin IT: Microarray in Ecological Research: A Case Study of a cDNA Microarray for Plant-Herbivore Interactions. BMS Ecology 2004., 4(13):
- Xie Y, Jeong KS, Pan W, Khodursky A, Carlin BP: A Case Study on Choosing Normalization Methods and Test Statistics for Two-Channel Microarray Data. Comparative and Functional Genomics 2004, 5: 432–444. 10.1002/cfg.416PubMed CentralView ArticlePubMedGoogle Scholar
- Bolstad B, Irizarry R, Astrand M, Speed T: A Comparison of Normalization Methods for High Density Oligonucleotide Array Data Based on Variance and Bias. Bioinformatics 2003, 19(2):185–193. 10.1093/bioinformatics/19.2.185View ArticlePubMedGoogle Scholar
- Dudoit S, Gentleman RC, Quackenbush J: Open Source Software for the Analysis of Microarray Data. BioTechniques 2003, 34: s45-s51.Google Scholar
- Saeed A, Sharov V, White J, Li J, Liang W, Bhagabati N, Braisted J, Klapa M, Currier T, Thiagarajan M, Sturn A, Snuffin M, Rezantsev A, Popov D, Ryltsov A, Kostukovich E, Borisovsky I, Liu Z, Vinsavich A, Trush V, Quackenbush J: TM4: A Free, Open-Source System for Microarray Data Management and Analysis. Bio Techniques 2003, 34: 374–378.Google Scholar
- Quackenbush J: Microarray Data Normalization and Transformation. Nature Genetics Supplement 2002, 32: 496–501. 10.1038/ng1032View ArticleGoogle Scholar
- Williams RD, King SN, Greer BT, Whiteford CC, Wei JS, Natrajan R, Kelsey A, Rogers S, Campbell C, Pritchard-Jones K, Khan J: Prognostic Classification of Relapsing Favorable Histology Wilms Tumor using cDNA Microarray Expression Profiling and Support Vector Machines. Genes, Chromosomes and Cancer 2004, 41: 65–79. 10.1002/gcc.20060View ArticlePubMedGoogle Scholar
- Zhu X, Hart R, Chang MS, Kim JW, Lee SY, Cao YA, Mock D, Ke E, Saunders B, Alexander A, Grossoehme J, Lin KM, Yan Z, Hsueh R, Lee J, Scheuermann RH, Fruman DA, Seaman W, Subramaniam S, Sternweis P, Simon MI, Choi S: Analysis of the Major Patterns of B Cell Gene Expression Changes in Response to Short-Term Stimulation with 33 Single Ligands. The Journal of Immunology 2004, 173: 7141–7149.View ArticlePubMedGoogle Scholar
- Khaitovich P, Weiss G, Lachmann M, Hellmann I, Enard W, Muetzel B, Wirkner U, Ansorge W, Paabo S: A Neutral Model of Transcriptome Evolution. PLoS Biology 2004, 2(5):682–689. 10.1371/journal.pbio.0020132View ArticleGoogle Scholar
- Watkinson JI, Sioson AA, Vasquez-Robinet C, Shukla M, Kumar D, Ellis M, Heath LS, Ramakrishnan N, Chevone B, Watson LT, van Zyl L, Egertsdotter U, Sederoff RR, Grene R: Photosynthetic Acclimation is Reflected in Specific Patterns of Gene Expression in Drought-Stressed Loblolly Pine. Plant Physiology 2003, 133(4):1702–1716. 10.1104/pp.103.026914PubMed CentralView ArticlePubMedGoogle Scholar
- Sioson AA, Watkinson JI, Vasquez-Robinet C, Ellis M, Shukla M, Kumar D, Ramakrishnan N, Heath LS, Grene R, Chevone BI, Kadafar K, Watson LT: Expresso and Chips: Creating a Next Generation Microarray Experiment Management System. In Proceedings of the Next Generation Software Systems Workshop, 17th International Parallel and Distributed Processing Symposium (IPDPS '03). Nice, France; 2003:209b.Google Scholar
- Heath LS, Ramakrishnan N, Sederoff RR, Whetten RW, Chevone BI, Struble CA, Jouenne VY, Chen D, van Zyl LM, Grene R: Studying the Functional Genomics of Stress Responses in Loblolly Pine using the Expresso Microarray Management System. Comparative and Functional Genomics 2002, 3: 226–243. 10.1002/cfg.169PubMed CentralView ArticlePubMedGoogle Scholar
- Galbraith D: Arabidopsis Oligonucleotide Microarrays.[http://www.ag.arizona.edu/microarray/]
- Mane SP, Vasquez-Robinet C, Sioson AA, Heath LS, Grene R: Phospholipase D alpha is Involved in Drought Stress Signaling in Arabidopsis . Poster presented at the International Conference on Plant Lipid-Mediated Signaling, Raleigh, NC 2005.Google Scholar
- Li P, Sioson AA, Mane SP, Ulanov A, Grothaus G, Heath LS, Murali TM, Bohnert HJ, Grene R: Response Diversity of Arabidopsis thaliana Ecotypes and Thellungiella halophila in elevated CO 2 in the field. Manuscript submitted 2005.Google Scholar
- Yang YH, Dudoit S, Luu P, Lin DM, Peng V, Ngai J, Speed TP: Normalization for cDNA Microarray Data: A Robust Composite Method Addressing Single and Multiple Slide Systematic Variation. Nucleic Acids Research 2002, 30(4):e15. 10.1093/nar/30.4.e15PubMed CentralView ArticlePubMedGoogle Scholar
- Yang I, Chen E, Hasseman J, Liang W, Frank B, Wang S, Sharov V, Saeed A, White J, Li J, Lee N, Yeatman T, Quackenbush J: Within the Fold: Assessing Differential Expression Measures and Reproducibility in Microarray Assays. Genome Biology 2002, 3(11):1–12.Google Scholar
- Cleveland W: Robust Locally Weighted Regression and Smoothing Scatterplots. J Amer Stat Assoc 1979, 74: 829–836. 10.2307/2286407View ArticleGoogle 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):S96-S104.View ArticlePubMedGoogle 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.