How to decide? Different methods of calculating gene expression from short oligonucleotide array data will give different results
© Millenaar et al; licensee BioMed Central Ltd. 2006
Received: 23 June 2005
Accepted: 15 March 2006
Published: 15 March 2006
Short oligonucleotide arrays for transcript profiling have been available for several years. Generally, raw data from these arrays are analysed with the aid of the Microarray Analysis Suite or GeneChip Operating Software (MAS or GCOS) from Affymetrix. Recently, more methods to analyse the raw data have become available. Ideally all these methods should come up with more or less the same results. We set out to evaluate the different methods and include work on our own data set, in order to test which method gives the most reliable results.
Calculating gene expression with 6 different algorithms (MAS5, dChip PMMM, dChip PM, RMA, GC-RMA and PDNN) using the same (Arabidopsis) data, results in different calculated gene expression levels. Consequently, depending on the method used, different genes will be identified as differentially regulated. Surprisingly, there was only 27 to 36% overlap between the different methods. Furthermore, 47.5% of the genes/probe sets showed good correlation between the mismatch and perfect match intensities.
After comparing six algorithms, RMA gave the most reproducible results and showed the highest correlation coefficients with Real Time RT-PCR data on genes identified as differentially expressed by all methods. However, we were not able to verify, by Real Time RT-PCR, the microarray results for most genes that were solely calculated by RMA. Furthermore, we conclude that subtraction of the mismatch intensity from the perfect match intensity results most likely in a significant underestimation for at least 47.5% of the expression values. Not one algorithm produced significant expression values for genes present in quantities below 1 pmol. If the only purpose of the microarray experiment is to find new candidate genes, and too many genes are found, then mutual exclusion of the genes predicted by contrasting methods can be used to narrow down the list of new candidate genes by 64 to 73%.
Affymetrix GeneChip arrays have become a widely used microarray platform. An advantage of such a broadly accepted platform is that it greatly helps to compare RNA expression patterns between experiments from different researchers. RNA expression patterns are a starting point for high-level analysis such as clustering, self-organising maps or metabolic pathway analysis. The importance of "low level analysis", that is, calculating the expression levels and the normalization procedure as explained below, is easily underestimated.
Calculating gene expression is more complex with Affymetrix arrays compared to spotted arrays, as the latest generation Affymetrix arrays have 11 probe pairs available to interrogate each gene, whereas spotted arrays in general have only one cDNA or oligo probe per gene. Ideally, the signal intensity from all of these 11 probes should be equal because they all interrogate the same gene. In reality, however, there are enormous differences between individual probes in a probe set. The pattern of the probe signals in a probe set is called a "probe response pattern" . The probe response pattern of one gene is roughly the same from all RNA samples isolated from a given tissue, irrespective of the treatment . The Affymetrix system uses Microarray Suite (MAS) or its successor GeneChip Operating Software (GCOS) for operation of the microarray scanner and also for calculation of the expression values. The expression values are calculated according to the so-called "One-step Tukey's biweight algorithm", which is based on a weighed median . Recently, more methods have become available for calculating gene expression values, such as "Model-Based Expression Indexes" (MBEI) from Li and Wong [1, 3] that are implemented in the freeware program dChip. Also available is the "Robust Multi-array Average" (RMA) software from Irizarry et al. [4, 5]. In contrast to MAS 5.0 or GCOS, in which information from only one microarray is used, model-based algorithms incorporate information from multiple microarrays to calculate the expression of a gene. The probe response pattern is fitted over multiple arrays with a multiplicative model in dChip and an additive model in RMA software. These algorithms use the fitted models to detect abnormally behaving probes, which are subsequently excluded for calculating gene expression. Therefore, gene expression from these model-based algorithms can be expected to provide more reproducible results. Both dChip and RMA use a stochastic model to estimate gene expression. A potentially better way of estimating gene expression would be the use of a physical model, that is, a model that explains the probe response pattern by physical laws. A few attempts have been made to make such physical models  such as "Positional-Dependent-Nearest-Neighbour model" (PDNN, ), or a combination of the two model types GeneChip and RMA (GC-RMA, ).
A specific feature of Affymetrix microarrays are the so-called MisMatch (MM) probes. Each of the 11 probe pairs in a probe set has a Perfect Match (PM) and a MisMatch probe, with each probe having 25 bases. The PM probes are designed to bind perfectly to the gene of interest and the MM probes have a contrasting base at position 13 with the intention to measure non-specific binding . The PM and MM probe of one probe pair are located physically next to each other on the microarray. Both MAS5 and dChip in the PMMM mode use this information to calculate gene expression. However, in the literature there is a debate over the usefulness of MM probes, which has resulted in new algorithms that use only the PM signals to calculate gene expression. These algorithms are implemented in dChip PM mode [1, 3], in the "Robust Multi-array Average" (RMA) software [4, 5], GC-RMA  and in PDNN from Zhang et al. .
When dealing with experiments involving multiple arrays, there is a need for data normalization, after calculating gene expression. The overall signal intensity can vary between arrays and if this variation does not have a biological origin, it should be removed. Non-biological variation can arise due to differences in sample preparation, production and processing of arrays. There are different algorithms available to normalize microarray data. The MAS 5.0 or GCOS software currently assumes that intensity differences between two or more arrays are linearly related and have a zero y-intercept. This allows for a very simple and robust normalization factor. However, the drawback is that this method cannot adjust for non-linear relations. Schadt et al.  found 10 to 50% differences in the slope between the low and high expressed genes in many arrays. In those cases an alternative normalization method is needed. For example, dChip uses a non-linear smoothing method . In MAS5 and dChip one array is chosen as the base line array, to which the other arrays are normalized. RMA does not use a base line array, but "quantile normalization", based on the complete data set . Furthermore, it is preferred that genes changing in expression should not be included in the normalization procedure. Therefore, some algorithms (like dChip) use invariant sets; these comprise a set of probes that most likely have the same expression between arrays. An invariant set is determined by ranking probe intensities such that when a probe has approximately the same ranking number in two arrays, then this probe is likely to have the same expression in both arrays [1, 3].
Thus different algorithms calculate gene expression and normalize the data in completely different ways. We hypothesize that the signal of gene expression, and therefore the identity of genes that are significantly up/down regulated, depends on the algorithm used. To test this hypothesis a microarray experiment was performed with Arabidopsis thaliana accession Columbia-0 and various calculation methods were compared. Plants were exposed for 3 hours to elevated ethylene levels (5 μl l-1) or were transferred to low-light levels for 3 hours (from 200 to 15 μmol m-2 s-1). Both conditions induce an identical phenotype; the upward bending of the leaves by 20 degrees (hyponastic growth; ). The treated and control petioles were harvested at the same time with three biological replicates. The dataset is publicly available . These data were analysed with six different methods: MAS5 (or GCOS), dChip PMMM, dChip PM-only, RMA, GC-RMA and PDNN. Some other studies have compared different calculation methods. Barash et al.  compared only three algorithms with a very limited amount of microarrays. Furthermore, Choe et al.  compared many options in the Bioconductor's Affy package, including the calculation of the expression values, however, again only three algorithms were compared. Seo et al.  compared six methods, four of which are also used in this paper.
Different algorithms are optimised for different experimental data sets and hence to different signal to noise ratios. Thus, the best method can depend on the specific experiment, so there might not be a perfect method suitable for all conditions . It is therefore, important to test for each experiment which method is the most appropriate.
Since different calculated expression patterns were expected between algorithms, various criteria were designed and tested to compare the different methods, and to select the best method. 1: We compared the performance of the six methods on the spike-in data set of Affymetrix . 2: Theoretically, we expected that model-based algorithms would produce more reproducible results because aberrant probes are efficiently removed. Therefore, we compared the distribution of the coefficient of variation for three replicated microarrays between the algorithms. 3: On the basis of data from the literature or previous experiments we expected that specific genes are affected by a treatment or mutation (biological relevance). 4: Furthermore, the problems with the use of MM signals as discussed in the literature are mainly statistical. Based on the following arguments we would like to discuss the use of MM values. i: First, we expected that cRNAs smaller than 13 base pairs bind perfectly to both probes. Calculation of the melting temperature is performed to estimate the likelihood of those events. ii: Secondly, we expected no correlation between the signals from the MM probes and the PM probes because the non-specific binding should have a random influence on the MM probes. To test this, we calculated correlation coefficients between PM and MM signals of 200 randomly chosen probe sets. 5: Microarray-obtained gene expression of a subset of the common genes was compared with Real-Time RT-PCR. Furthermore, we compared the expression from genes which are exclusively predicted with the best method.
The proposed criteria for biological relevance and reproducibility are easily applied by any microarray user on their own dataset. On the basis of all of the criteria described above a choice for the most suitable algorithm for each dataset is rationalized.
Results and discussion
Differences between six array algorithms
We hypothesized that the calculated level of gene expression is dependent on the algorithm used because these are based on completely different ways of calculating gene expression and normalization. As a result, a different set of genes may be identified as differentially expressed between a treatment and control.
The number of differentially expressed genes. Number of genes significant (Ttest p < 0.05) up or down regulated between air control and ethylene or low-light treatment in six different algorithms. The number of biological replicates for each treatment is three. The number of probe pairs per chip is 22746.
Control – ethylene
Control – Low-light
When instead of a T-test with p < 0.05 more stringent conditions were used (p < 0.01; < 0.001; < 0.0001), less genes were found to be significantly regulated (2483; 798; 137; 21, respectively) and also the percentage genes significantly regulated in all algorithms decreased (32%; 16%; 2%; 0%). Moreover, if genes were selected that showed more than twofold up or down regulation, different numbers were obtained. However, the huge differences between the algorithms remained (data not shown). The set of genes that was detected by all methods had on average a lower p value (0.013) than the genes identified solely by one method (p = 0.028). On average the percentage genes that overlapped was 61% over all possible double combinations but dropped to 48 (all triple combinations), 40 (4), 35 (5) and finally to 32% (6). Our results showed a poor overlap between algorithms. This is not in full agreement with Barash et al.  and Choe et al.  who showed 60 to 70% overlap between algorithms. However, their results were based on a small number of microarrays or on a special (only) spike-in experiment. To our knowledge this is the first time that such large differences between algorithms are reported.
From the Venn diagram (Fig. 1) it is clear that MAS5 has the highest absolute number of uniquely identified differentially expressed genes compared to the other methods. When expressed relatively, MAS5 had the highest percentage of unique genes (26%), followed by PMMM (18%), PDNN (12%), PM (12%), RMA (10%) and GCRMA (9%). Furthermore, there are six possible combinations in which five methods might identify the same genes (Fig. 1; abcde, abcdf, abcef, abdef, acdef, bcdef), and the only combination in which MAS5 was not included (abcef) contains dramatically more genes (291) than the other five combinations in which MAS5 was included. These data show that MAS5 gives more unique results compared to the other five methods (Fig. 1).
Correlations between the expression values. Pearson correlation between the expression values calculated with six different methods. The average signal intensity from all 22746 probe sets from the three control chips was used.
There were large differences between the six algorithms with respect to both gene expression and the genes that were significantly up/down regulated. This was a consequence of different ways of calculating gene expression and normalization of the microarrays. The discrepancy between the different algorithms can also be exploited. When the only aim of a microarray experiment is to select new candidate genes then the list of new candidate genes is shortened by approximately 68% by using genes predicted by all six methods. However, when the aim is to obtain expression values from a microarray experiment, a decision regarding which method to use has to be made. Various criteria were designed and tested to compare the different methods. Cope et al.  compared different methods and focused their research on the assessments of performance in terms of bias and variance. This was studied on the spike-in experiment from Affymetrix with the HGU95A chips and on a RNA dilution study by GeneLogic. These authors also wrote a program to compare other available algorithms . Our study further explored the performance of individual algorithms further by five criteria: 1) comparison with spiked-in genes, 2) reproducibility, 3) biological relevance, 4) the use of MM probes, and 5) comparison with Real Time RT-PCR.
1) Spike-in comparison
A commonly used RNA spike-in experiment from Affymetrix was used to test which method produced the most accurate results. In this spike-in experiment RNA of 42 genes in 14 concentrations were added to a human RNA sample. Virtually all signal intensities found in real microarray experiments are covered by the signal intensities found for the spike-in genes. In this experiment, the relation between known concentrations of spiked-in RNA and signal intensity can be compared between the six methods.
As a measure of precision, the significance level was calculated between successive concentrations of all of the 42 RNA's spiked-in, so between 0 and 0.125 pM, between 0.125 and 0.25 pM up to 512 pM. Subsequently the fraction of genes that were up regulated significantly was plotted (figure 4B). For the higher concentrations of spike-in genes (>1 pmol RNA or zero on a ln scale), dChip PM, dChip PMMM, RMA and GC-RMA showed almost perfect results, this in contrast to PDNN and MAS5. The results of the ROC curves from Cope et al.  also showed that MAS5 performed worse than RMA and PM. For the genes expressed at low levels, all methods performed inadequately. Also, when all successive spike-in concentrations were compared with the zero control, again the six methods showed adequate results only at higher concentrations (>1 pmol, data not shown). This shows that although there are differences between methods in accuracy at low concentrations (Fig. 4A), these do not seem to be relevant, since no significant differences between any of the methods were observed in this lower range (Fig. 4B).
Minimum signal. Signal intensity from the RNA spike-in experiments at a concentration (1 pmol) where the spiked-in concentrations were significant different from the background. The number (#) and percentage (%) genes below the minimal signal (Min. signal) is obtained from non-transformed data.
3) Biological relevance
Biological comparisons. Expression of genes in ethylene biosynthesis and signal transduction cascade as calculated with six different methods. Numbers in brackets after a gene name show how many genes of a family are on the microarray, followed by the size of the family. Expression data of a gene family is the summation of expression of the individual gene members. The columns with a p showing the significance of the fold change ('x' columns) between air control and ethylene (eth) treated plants. ns = not significant * = p < 0.05, ** = p < 0.01, *** = p < 0.001.
ACC synthase (10/11)
ACC oxidase (10/10)
ACC synthase (10/11)
ACC oxidase (10/10)
4) The usefulness of mismatch (MM) probes
One possible contributor to the MM values may be binding of small cRNAs to the lateral sequence of both PM or MM probes. In other words, since the difference between the PM and MM is at position 13 of an oligonucleotide 25 bases long, any cRNA that can bind to either of those lateral 12 bases might remain bound in the detection step. It is similarly likely that a cRNA fragment with its terminal 12 base pairs complementary to either lateral 12 base sequence of the probe may also remain bound at detection. In this latter case if the cRNA represents a truncated version of the expected signal cRNA then true transcript signal may be added to both PM and MM values. According to the rules regarding maximum GC content for probes described by Lockhart et al.  the experimental hybridisation for affymetrix arrays is performed at a temperature (45°C) and stringency wash (50°C) which are both higher than the maximum possible melting temperature of a 12 mer. Based on this argument, cRNAs smaller than 13 bases will not contribute to the microarray signal values. Furthermore, since the cRNAs are fragmented before hybridisation to an average length of 100 bases long the percentage of small stretches is low.
If MM signals are important for detecting genes expressed at low levels, differences are expected between algorithms using the MM signals (MAS, PMMM) and methods not using the MM signal. Analysis of the spike-in data shows (Fig. 4B) that MAS and PMMM are, at best, performing equally well in detecting low expressed spiked-in genes. Of the four methods that do not use the MM signals one performed poorly for highly expressed spike-in genes, compared to one out of two for the methods using the MM signals. These data suggest that the MM signals do not increase the detection of spike-in genes in this data set. Future research with independent techniques and independent data sets will have to prove how useful MM signals are.
5) Real Time RT-PCR
To further test which algorithm calculates microarray gene expression most accurately, the results of the six algorithms were compared with Real Time RT-PCR on a subset of 18 genes. These 18 genes are predicted by all the six methods to be differentially regulated. Real Time RT-PCR data can also be calculated with several different algorithms. Most widely used is the ΔΔCt method . One of the assumptions of this method is that the efficiencies of the PCR reactions are the same for the gene of interest and for a reference gene. More recently, other methods have become available, such as the "Assumption-free analysis" . This method determines the efficiency of each PCR reaction by calculating the slope for the exponential part of the curve. Czechowski et al.  combined the two methods by using the average efficiency of a gene to correct each ΔΔCt value.
Correlation between microarray and RT-PCR data. Pearson correlation of expression values between six microarray algorithms and three Real Time RT-PCR algorithms on a subset of 18 genes. The negative sign from the numbers in the second and last column are removed, because a lower ΔCt value means more gene product present, consequently leading to a negative correlation with microarray data. All correlations were significant. ns = not significant (p < 0.01). * = p < 0.05, ** = p < 0.01
ΔCt + Assumption Free
To test if there are differences between genes (18) with respect to what the best method for calculating expression levels, the same correlations were calculated per gene. However, due to the low number of data points per correlation (3), hardly any significant correlations were found. Therefore, we were not able to determine if there is a gene dependency for the best method. To resolve this matter, a wider range of expression levels per gene should be examined.
We also tested if the correlations between the microarray and Real Time RT-PCR data depended on signal intensity. Possibly such correlations might be weak or absent at low signal intensities. However, we did not find any significant relation between these parameters. Moreover, except for one, the signal intensities where all larger than the minimal signal intensity to detect spike in gene expression as found in table 3.
RT-PCR data of RMA specific genes. Ethylene/air expression signals from 10 genes calculated with RMA and compared to the Real Time RT-PCR data as calculated with ΔCt. All RMA signals from ethylene treated plants are significantly (0.05 > p > 0.01) different from the air controls (average p = 0.037). The Real Time RT-PCR data shown, are averages, standard deviations and p values after comparison with the control. The average control signal is set to 1. The numbers in bold indicate genes which show a change in expression similar for both methods. For RMA, n = 3 and for Real Time RT-PCR, n = 3–5.
Gene expression was calculated with six algorithms, MAS5, dChip PMMM, dChip PM, RMA, GC-RMA and PDNN. All six algorithms resulted in different levels of gene expression, with MAS5 calculating the lowest and PDNN the highest signals. Consequently, depending on the method used, different genes will be identified as up or down regulated by a given treatment. Therefore, it is up to the user to decide which method to use to calculate the gene expression profile. The spike-in experiment showed that the differences between the methods for the low expressed genes are not relevant since measurements up to 1 pmol RNA, did not result in significant signals above background. Furthermore, for relatively high expressed genes MAS5 and PDNN showed the lowest percentage of correctly calculated genes. In our comparison all six methods yielded the expected biological results, with respect to ethylene-induced gene expression. Furthermore, all model-based algorithms produced more reproducible results in comparison to MAS5. Moreover, we conclude that MM signals do not represent non-specific binding for PM signals, since in 47.5% of the cases there was a significant correlation between the MM and PM signals. Subtracting the MM from the PM signals will therefore most likely underestimate the expression signal. RMA generally yields the most reproducible results. However, we were not able to verify most genes which were solely identified as having changed expression by RMA with Real Time RT-PCR. Different experimental systems may require different methods to obtain the best result . If the only purpose of the microarray experiment is to find new candidate genes, and too many genes are found, then mutual exclusion of the genes predicted by contrasting methods can be used to narrow down that list of candidate genes.
Plant material and growth conditions
Arabidopsis thaliana accession Columbia-0 (Col-0; N1092, Nottingham Arabidopsis Stock Centre) was used. Seeds were sown on moistened filter paper in sealed Petri-dishes and cold-stratified in the dark at 4°C for 4 d. Subsequently, the seeds were germinated for 4 d in a growth chamber with the following conditions: 20°C, 70% (v/v) relative humidity, 9 h photoperiod: 200 μmol m-2 s-1 photosynthetic active radiation (PAR). Seedlings were transferred with a brush to pots (70 ml) containing a mixture of potting soil and perlite (1:2, v:v), enriched with 0.14 mg MgOCaO (17%; Vitasol BV, Stolwijk, The Netherlands) and 0.14 mg slow-release fertilizer (Osmocote "plus mini"; Scotts Europe bv, Heerlen, The Netherlands) per pot. Prior to seedling transfer, to each pot 20 ml nutrient solution was added, containing: 2.6 mM KNO3, 2.0 mM Ca [NO3]2, 0.6 mM KH2PO4, 0.9 mM MgSO4, 6.6 μM MnSO4, 2.8 μM ZnSO4, 0.5 μM CuSO4, 66 μM H3BO3, 0.8 μM Na2MoO4 and 134 μM Fe-EDTA, pH 5.8. All chemicals were p.a. grade and obtained from Merck (Darmstadt, Germany). Following transplantation, plants were grown for 28 d in a growth chamber (conditions as described above). Pots with seedlings were kept in a glass-covered tray for the first 4 d following transplantation, after which they were transferred to irrigation mats (Maasmond-Westland, The Netherlands). The mats were automatically watered with tap water to saturation once a day (at the beginning of the light period), and the excess water was drained.
Ethylene and low-light experiments
Ethylene (100 μL L-1; Hoek Loos BV, The Netherlands) and air (70% relative humidity) were mixed using flow meters (Brooks Instruments BV, The Netherlands) to generate a concentration of 5 μL L-1 ethylene, which was flushed continuously through glass cuvettes (13.5 × 16.0 × 29.0 cm) at 75 L h-1, and then vented to the outside of the building. This concentration saturates ethylene-induced petiole elongation in R. palustris . A concentration of 1 μL L-1 ethylene was reached in the cuvettes after approximately 10 min of starting the treatment; 5 μL L-1 was reached after 40 min. The ethylene concentration was checked regularly on a gas chromatograph (GC955, Synspec, Groningen, Netherlands), and remained constant for the duration of the experiment. Control cuvettes were flushed with air (70% relative humidity) at the same flow rate.
For the shading treatment, the light quantity was reduced by 90% to 15 – 20 μmol m-2 s-1 (photosynthetic active radiation) at the start of the experiment. This was achieved by switching off a number of lamps in the growth chamber and by using a mesh black spectrally neutral shade cloth. These alterations did not change the light quality when checked with a Licor1800 spectro-radiometer (Li-cor, Lincoln, Ne, USA) (data not shown).
Petioles were harvested and directly frozen in liquid nitrogen from plants in stage 3.9 according to Boyes et al. , which means a plant with approximate 17 rosette leaves. Subsequently RNA was isolated from these petioles with RNeasy extraction protocol from Qiagen (Valencia, CA, USA).
Approximately 10 μg of total RNA was reverse transcribed at 42°C for 1 h to generate first strand DNA using 100 pmol oligo-dT(24) primer containing a 5'-T7 RNA polymerase promoter sequence (5'-GGCCAGTGAATTGTAATACGACTCACTATAGGGAGGCGG-(dT)24-3'), 50 mM Tris-HCl (pH 8.3), 75 mM KCl, 3 mM MgCl2, 10 mM dithiothreitol (DTT), 10 mM dNTPs and 200 units SuperScript II reverse transcriptase (Invitrogen Life Technologies, Breda, The Netherlands). Following first strand synthesis, second strand DNA was synthesized using 10 units of E. coli polymerase I, 10 units of E. coli DNA ligase and 2 units of RNase H in a reaction containing 25 mM Tris-HCl (pH 7.5), 100 mM KCl, 5 mM MgCl2, 10 mM (NH4)SO4, 0.15 mM b-NAD+ and 10 mM dNTPs. The second strand synthesis reaction proceeded at 16°C for 2 hours before 10 units of T4 DNA polymerase was added and the reaction allowed to proceeded for 5 minutes. The reaction was terminated by adding 0.5 M EDTA. Double stranded cDNA products were purified using the GeneChip Sample Cleanup Module (Affymetrix, USA).
The synthesized cDNAs were in vitro transcribed by T7 RNA polymerase (ENZO BioArray High Yield RNA Transcript Labeling Kit, Enzo, San Francisco, CA, USA) using biotinylated nucleotides to generated complementary RNAs (cRNAs). The cRNAs were purified using the GeneChip Sample Cleanup Module (Affymetrix, USA). The cRNAs were then randomly fragmented at 94°C for 35 min in a buffer containing 40 mM Tris-acetate (pH 8.1), 100 mM potassium acetate, and 30 mM magnesium acetate to generate molecules of approximately 35 to 200 bases long.
The fragmented cRNA were mixed with 0.1 mg ml-1 of sonicated herring sperm DNA in a hybridisation buffer containing 100 mM 2-N-morpholino-ethane-sulfonic acid (MES), 1 M NaCl, 20 mM EDTA and 10% Tween-20 making a hybridisation mixture. The hybridisation mixture containing the fragmented cRNA was denatured at 99°C for 5 min and equilibrated at 45°C for 5 min and centrifuged at maximum speed in a microcentrifuge for 5 min to remove any insoluble material from the hybridisation mixture. The hybridisation mix was transferred to the Arabidopsis ATH1-121501 genome array (Affymetrix) cartridge and hybridised at 45°C for 16 h on a rotisserie at 60 rpm. After a 16 h hybridisation period, the arrays were washed and stained in a fluidics station (Affymetrix). The arrays were initially washed in a low stringency buffer A (6× SSPE (0.9 m NaCl, 0.06 M NaH2PO4, 0.006 M EDTA), 10% Tween 20) at 25°C for 10 min and then incubated with a high stringency buffer B (100 mM MES, 0.1 NaCl, 10% Tween-20) at 50°C for 20 min, then stained with 10 mg ml-1 of Streptavidin Phycoerythrin (SAPE), in stain buffer containing 100 mM MES, 1 M NaCl, 0.05% Tween 20 and 2 mg ml-1 BSA at 25°C for 10 min, washed with wash buffer A at 25°C for 20 min and stained with biotinylated anti-streptavidin antibody at 25°C for 10 min. After antibody staining the arrays were stained again with SAPE for signal amplification and washed with buffer A at 30°C for 30 min. The arrays were finally scanned and the intensities averaged with a Agilent GeneArray Scanner.
Real time RT-PCR
For one RNA sample eight petioles obtained from two plants were mixed and ground in liquid nitrogen. Samples were taken after 3 h of ethylene treatment or from plants in air for 3 h. Total RNA was isolated from Arabidopsis thaliana petioles using RNeasy Plant Mini Kit (Valencia, CA, USA). Genomic DNA was removed using the DNA-Free kit (Ambion, Cambridgeshire, United Kingdom). cDNA was synthesized using 3.3 μg total RNA with Superscript III RNase H- Reverse Transcriptase (Invitrogen, Breda, The Netherlands) using Random-Hexamer Primers. Real time PCR reactions were performed on MyiQ Single-Color Real-Time PCR Detection System and Software using iQ SYBR Green Supermix fluorescein (Bio-Rad laboratories, Veenendaal, The Netherlands).
RT-PCR primers. Primers for the genes used in Real-Time RT-PCR and the length of the PCR product. Each gene is represented with the Affymetrix probe number and the Arabidopsis Genome Initiative (AGI) code.
Genes common for all microarray methods
RMA specific genes
The RNA spike-in experiment is carried out by Affymetrix on Human HG-U133 chips. The CEL files are available online . In this experiment the RNA of 42 genes are spiked-in at 14 concentrations 0, 0.125, 0.25, 0.5, 1, 2, 4, 8, 16, 32, 64, 128, 256, 512 pmol and in total 42 micro arrays were used. For more details see the Affymetrix website.
The following programs are used, MAS 5.0 from Affymetrix (Santa Clara, USA), DNA-Chip analyser (dChip) version 1.2 , RMA version 0.2 , PDNN version 2.3 , GC-RMA as was incorporated in the program ArrayAssist Lite version 3.2.1954.31591 (Stratagene, La Jolla, CA, USA) and Excel 2000 (Microsoft, Redmond, USA). Cel files from MAS5 containing the signal intensity per feature were used as input for dChip, RMA and PDNN. A two-sided T-test with unequal variances was performed in Excel with a p value of 0.05 unless stated otherwise. The coefficient of variation or CV is calculated in Excel as the standard deviation divided by the average signal. The CV was sorted separately for each algorithm before the distribution of the CV was plotted; this was also the procedure for the signal distributions, slope and correlation data in figure 2, 5, and 7. A macro was written in Excel to find common genes between two sets of genes to construct figure 1, this macro is available .
Real Time RT-PCR data were calculated in three ways: 1) comparative Ct method described by Livak and Schmittgen , 2) "Assumption-free analysis" by Rademakers et al.  and 3) a combination of the first two methods by Czechowski et al. .
This work was supported by a PIONIER grant (800.84.470) of the Dutch Science Foundation (NWO) to LACJV and by a NIELS STENSEN grant to FFM. Furthermore, we thank Dr. Ronald Pierik, Dr. Frank C.P. Holstege (Utrecht University, The Netherlands) and Dr. Alex Boonman (University of Amsterdam, The Netherlands) for their comments on earlier versions of this manuscript, and Dr. Matthijs J. Coster for helping to design the Venn diagram.
- Li C, Wong WH: Model-based analysis of oligonucleotide arrays: expression index computation and outlier detection. Proc Nat Acad Sci USA 2001, 98: 31–36. 10.1073/pnas.011404098PubMed CentralView ArticlePubMedGoogle Scholar
- Affymetrix: Microarray Suite User Guide. Affymetrix 2001., Version 5: [http://www.affymetrix.com/support/technical/manuals.affx]Google Scholar
- Li C, Wong WH: Model-based analysis of oligonucleotide arrays: model validation, design issues and standard error application. Genome Biology 2001, 8: 0032.Google Scholar
- Irizarry RA, Bolstad BM, Collin F, Cope LM, Hobbs B, Speed TP: Summaries of Affymetrix genechip probe level data. Nucleic Acids Research 2003, 31: e15. 10.1093/nar/gng015PubMed CentralView ArticlePubMedGoogle Scholar
- Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP: Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics 2003, 4: 249–264. 10.1093/biostatistics/4.2.249View ArticlePubMedGoogle Scholar
- Hekstra D, Taussig AR, Magnasco M, Naef F: Absolute mRNA concentrations from sequence-specific calibration of oligonucleotide arrays. Nucleic Acids Research 2003, 31: 1962–1968. 10.1093/nar/gkg283PubMed CentralView ArticlePubMedGoogle Scholar
- Zhang L, Miles MF, Aldape KD: A model of molecular interactions on short oligonucleotide microarrays. Nature Biotechnology 2003, 21: 818–821. 10.1038/nbt836View ArticlePubMedGoogle Scholar
- Wu Z, Irizarry RA: Preprocessing of oligonucleotide array data. Nature Biotechnology 2004, 22: 656–658. 10.1038/nbt0604-656bView ArticlePubMedGoogle Scholar
- Lipshutz RJ, Fodor SPA, Gingeras TR, Lockhart DJ: High density synthetic oligonucleotide arrays. Nature Genetics 1999, (Suppl 21):20–24. 10.1038/4447Google Scholar
- Schadt EC, Li C, Ellis B, Wong WH: Feature extraction and normalization algorithms for high-density oligonucleotide gene expression array data. J Cell Biochem 2001, (Suppl 37):120–125. 10.1002/jcb.10073Google 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 2003, 19: 185–193. 10.1093/bioinformatics/19.2.185View ArticlePubMedGoogle Scholar
- Millenaar FF, Cox MCH, de Jong van Berkel YEM, Welschen RAM, Pierik R, Voesenek LAJC, Peeters AJM: Ethylene-induced differential growth of petioles in arabidopsis. Analyzing natural variation, response kinectics, and regulation. Plant Physiology 2005, 137: 998–1008. 10.1104/pp.104.053967PubMed CentralView ArticlePubMedGoogle Scholar
- Barasch Y, Dehan E, Krupsky M, Franklin W, Geraci M, Friedman N, Kaminski N: Comparative analysis of algorithms for signal quantitation from oligonucleotide microarrays. Bioinformatics 2004, 20: 839–846. 10.1093/bioinformatics/btg487View ArticleGoogle Scholar
- Choe SE, Boutros M, Michelson AM, Church GM, Halfon MS: Preferred analysis methods for Affymetrix GeneChips revealed by a wholly defined control dataset. Genome Biology 2005, 6: R16. 10.1186/gb-2005-6-2-r16PubMed CentralView ArticlePubMedGoogle Scholar
- Seo J, Bakay M, Chen Y-W, Hilmer S, Sheiderman B, Hoffman EP: Interactively optimizing signal-to-noise ratios in expression profiling: project-specific algorithm selection and detection p-value weighting in Affymetrix microarrays. Bioinformatics 2004, 20: 2534–2544. 10.1093/bioinformatics/bth280View ArticlePubMedGoogle Scholar
- Cope LM, Irizarry RA, Jaffee HA, Wu Z, Speed TP: A benchmark for Affymetrix genechip expression measures. Bioinformatics 2004, 20: 323–331. 10.1093/bioinformatics/btg410View ArticlePubMedGoogle Scholar
- Kieber JJ, Rothenberg M, Roman G, Feldmann KA, Ecker JR: CTR1, a negative regulator of the ethylene response pathway in Arabidopsis, encodes a member of the raf family of protein kinases. Cell 1993, 72: 427–441. 10.1016/0092-8674(93)90119-BView ArticlePubMedGoogle Scholar
- Peck SC, Pawlowski K, Kende H: Asymmetric responsiveness to ethylene mediates cell elongation in the apical hook of peas. The Plant Cell 1998, 10: 713–719. 10.1105/tpc.10.5.713PubMed CentralView ArticleGoogle Scholar
- Vriezen WH, Hulzink R, Mariani C, Voesenek LACJ: 1-Aminocyclopropane-1-carboxylate oxidase activity limits ethylene biosynthesis in Rumex palustris during submergence. Plant Physiology 1999, 121: 189–195. 10.1104/pp.121.1.189PubMed CentralView ArticlePubMedGoogle Scholar
- Wilkinson JQ, Lanahan MB, Yen H-C, Giovannoni JJ, Klee HJ: An ethylene-inducible component of signal transduction encoded by Never-ripe. Science 1995, 270: 1807–1809.View ArticlePubMedGoogle Scholar
- Vriezen WH, Van Rijn CPE, Voesenek LACJ, Mariani C: A homologue of the Arabidopsis thaliana ERS gene is actively regulated in Rumex palustris upon flooding. Plant journal 1997, 11: 1265–1271. 10.1046/j.1365-313X.1997.11061265.xView ArticlePubMedGoogle Scholar
- Hua J, Sakai H, Nourizadeh S, Chen QG, Bleecker AB, Ecker JR, Meyerowitz EM: EIN4 and ERS2 are members of the putative ethylene receptor gene family in Arabidopsis. Plant Cell 1998, 10: 1321–1332. 10.1105/tpc.10.8.1321PubMed CentralView ArticlePubMedGoogle Scholar
- Leclercq J, Adams-Phillips LC, Zegzouti H, Jones B, Latché A, Giovannoni JJ, Pech J-C, Bouzayen M: LeCTR1, a Tomato CTR1-Like Gene, Demonstrates Ethylene Signaling Ability in Arabidopsis and Novel Expression Patterns in Tomato. Plant Physiology 2002, 130: 1132–1142. 10.1104/pp.009415PubMed CentralView ArticlePubMedGoogle Scholar
- Lee I, Dombkowski AA, Athey BD: Guidelines for incorporating non-perfectly matched oligonucleotides into target-specific hybridization probes for a DNA microarray. Nucleic Acids Research 2004, 32: 681–690. 10.1093/nar/gkh196PubMed CentralView ArticlePubMedGoogle Scholar
- Naef F, Hacker CR, Patil N, Magnasco M: Empirical characterization of the expression ratios noise structure in high-density oligonucleotide arrays. Genome Biology 2002, 3: 0018. 10.1186/gb-2002-3-4-research0018View ArticleGoogle Scholar
- Naef F, Lim DA, Patil N, Magnasco M: DNA hydridization to mismatched templates: A chip study. Physical Review E 2002, 65: 040902. 10.1103/PhysRevE.65.040902View ArticleGoogle Scholar
- Zhou Y, Abagyan R: Match-only integral distribution (MOID) algorithm for high-density oligonucleotide array analysis. BMC Bioinformatics 2002, 3: 3. 10.1186/1471-2105-3-3PubMed CentralView ArticlePubMedGoogle Scholar
- Allemeersch J, Durinck S, Vanderhaeghen R, Alard P, Maes R, Seeuws K, Bogaert T, Coddens K, Deschouwer K, Van Hummelen P, Vuylsteke M, Moreau Y, Kwekkeboom J, Wijfjes AHM, May S, Beynon J, Hilson P, Kuiper MTR: Benchmarking the CATMA microarray. A novel tool for Arabidopsis transcriptome analysis. Plant Physiology 2005, 137: 588–601. 10.1104/pp.104.051300PubMed CentralView ArticlePubMedGoogle Scholar
- Lockhart DJ, Dong H, Byrne MC, Follettie MT, Gallo MV, Chee MS, Mittmann M, Wang C, Kobayashi M, Horton H, Brown EL: Expression monitoring by hydridization to high-density oligonucleotide arrays. Nature Biotechnology 1996, 14: 1675–1680. 10.1038/nbt1296-1675View ArticlePubMedGoogle Scholar
- Chee M, Yang R, Hubbell E, Berno A, Huang XC, Stern D, Winkler J, Lockhart DJ, Morris MS, Fodor SPA: Accessing Genetic information with high-density DNA arrays. Science 1996, 274: 610–614. 10.1126/science.274.5287.610View ArticlePubMedGoogle Scholar
- Wang DG, Fan J-B, Siao C-J, Berno A, Young P, Sapolsky R, Ghandour G, Perkins N, Winchester E, Spencer J, Kruglyak L, Stein L, Hsie L, Topaloglou T, Hubbell E, Robinson E, Mittmann M, Morris MS, Shen N, Kilburn D, Rioux J, Nusbaum C, Rozen S, Hudson TJ, Lipshutz R, Chee M, Lander ES: Large-scale identification, mapping, and genotyping of single-nucleotide polymorphisms in the human genome. Science 1998, 280: 1077–1082. 10.1126/science.280.5366.1077View ArticlePubMedGoogle Scholar
- Livak KJ, Schmittgen TD: Analysis of relative gene expression data using real-time quantitative PCR and the 2 -ΔΔCT Method. Methods 2001, 25: 402–408. 10.1006/meth.2001.1262View ArticlePubMedGoogle Scholar
- Rademakers C, Ruijter JM, Lekanne Deprez RH, Moorman AFM: Assumption-free analysis of quantitative real-time polymerase chain reaction (PCR) data. Neuroscience Letters 2003, 339: 62–66. 10.1016/S0304-3940(02)01423-4View ArticleGoogle Scholar
- Czechowski T, Bari RP, Stitt M, Scheible W-R, Udvardi MK: Real-time RT-PCR profiling of over 1400 Arabidopsis transcription factors: unprecedented sensitivity reveals novel root- and shoot-specific genes. Plant Journal 2004, 38: 366–379. 10.1111/j.1365-313X.2004.02051.xView ArticlePubMedGoogle Scholar
- Sokal RR, Rohlf FJ: Biometry, the principles and practice of statistics in biological research. New York: W.H. Freeman and company; 1995:582.Google Scholar
- Voesenek LACJ, Blom CWPM: Growth responses of Rumex species in relation to submergence and ethylene. Plant Cell Environ 1989, 12: 433–439.View ArticleGoogle Scholar
- Boyes DC, Zayed AM, Ascenzi R, McCaskill AJ, Hoffman NE, Davis KR, Görlach J: Growth stage-based phenotypic analysis of Arabidopsis: a model for high throughput functional genomics in plants. Plant cell 2001, 13: 1499–1510. 10.1105/tpc.13.7.1499PubMed CentralView ArticlePubMedGoogle Scholar
- http://biosun1.harvard.edu/~cli/ or http://biosun1.harvard.edu/complab/dchip/
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.