Expression profiling is a universal tool, with a range of applications that benefit from the accurate determination of differential gene expression. To allow normalization using endogenous transcript levels, current microarray analyses assume that relatively few transcripts vary, or that any changes that occur are balanced. When normalization using endogenous genes is carried out, changes in expression levels are calculated relative to the behaviour of most of the transcripts. This does not reflect absolute changes if global shifts in messenger RNA populations occur. Using external RNA controls, we have set up microarray experiments to monitor global changes. The levels of most mRNAs were found to change during yeast stationary phase and human heat shock when external controls were included. Even small global changes had a significant effect on the number of genes reported as being differentially expressed. This suggests that global mRNA changes occur more frequently than is assumed at present, and shows that monitoring such effects may be important for the accurate determination of changes in gene expression.
Microarray expression profiling is a universal tool, with a range of applications that benefit from the accurate determination of differential gene expression (Brown & Botstein, 1999; Young, 2000). The detection of changes in messenger RNA expression requires normalization between samples. This counters non‐biological variation, such as differences in labelled material, local array differences, dye‐specific biases and so on (Quackenbush, 2001; Tseng et al., 2001; Yang et al., 2002).
Normalization methods consist of an algorithm and the features on the array to which the algorithm is applied. Ideally, such features should have identical signals in all the samples under investigation. Most researchers carrying out microarray analyses have dismissed the idea of using invariant house‐keeping genes that are stably expressed across a wide range of experimental conditions (Lee et al., 2002). With a few exceptions (Talaat et al., 2002; Yang et al., 2002), most microarray experiments make use of the expression levels of all genes as normalization features. The assumption underlying this ‘all‐genes’ approach is that relatively few transcript levels vary between samples, or that any changes that occur are balanced. Normalization using the expression levels of endogenous genes means that changes are calculated relative to the majority of transcripts. These relative changes do not reflect the absolute changes at the cellular level that occur if global shifts in mRNA populations take place.
The aim of this study was to set up microarray experiments, incorporating external controls, to monitor the effects of inactivating components of the generally required transcription machinery, such as RNA polymerase II (Holstege et al., 1998; Wang et al., 2002). The purpose of external normalization controls is to derive a set of signals for which the final outcome is known to be equal among samples. This can be achieved by the addition of equal amounts of control RNA molecules to samples before processing. Here, we present results from experiments in which we observed global changes occurring under conditions that, in the past, have been studied without the use of external controls. As well as showing how global effects can be monitored, the results suggest that global mRNA changes occur more frequently than is presently assumed by researchers carrying out microarray analyses. This has important implications for the interpretation of such experiments.
External normalization controls
A generalized scheme for carrying out microarray analyses using external controls for normalization is shown in Fig. 1A. In this study, total RNA was spiked with a mixture of nine control RNAs. The concentrations of the control RNAs were varied over three orders of magnitude to cover a range of mRNA expression levels (Fig. 1B). Oligonucleotide probes (of 70 nucleotides in length) representing each control were spotted at least twice onto the microarray subgrids to generate sufficient data points (960 in total; 20 per subgrid). This allows local, expression‐dependent or intensity‐dependent normalization (Yang et al., 2002). The microarrays also incorporated other controls in addition to the gene probes (6,371 for Saccharomyces cerevisiae, each spotted twice, and 16,735 for the human arrays).
Global changes during mammalian heat shock
One study that has shown global changes in the mRNA population is that of the heat‐shock response of primary human umbilical vein endothelial cells (HUVECs). This is similar to many previous microarray studies that have examined cellular responses, and was investigated as part of our programme to understand transcriptional regulation. When external controls were used, a change in global mRNA levels during heat shock was shown, as indicated by the gradual separation of the green gene spots away from the blue control spots (Fig. 2A–D). This culminated in an almost twofold median drop in mRNA levels on average in the dye‐swap experiments (Fig. 2D). Typically, these responses have previously been normalized using expression levels of all genes. Approximately equal numbers of genes were interpreted as being up‐ or downregulated after normalization of the experiment in this way (Fig. 2E, left column), which is a consequence of making the assumption that there is no overall change. Such changes in individual transcript levels are relative to the behaviour of most of the transcripts. If a global change occurs, as is the case here (Fig. 2D), normalization using genes does not correlate with changes at the cellular level. The interpretation of these experiments differs markedly when external controls are used (Fig. 2E, middle column). The number of genes reported as being downregulated increased from 506 to 6,872 when an arbitrary twofold threshold was used. These changes are relative to equivalent amounts of total RNA because spiking of the controls was carried out at the total RNA stage. Counting cells before harvesting allows differences in RNA content in the two states to be taken into account (Fig. 2E, right column).
Results of small global changes during serum starvation
A second example shows the results of subjecting human mammary gland adenocarcinoma (MCF7) cells to serum deprivation for 30 h (Fig. 3). Although this resulted in only a small change in global mRNA levels (1.3‐fold median drop; Fig. 3A), it had a significant effect on the outcome. More than twice as many genes were found to be differentially expressed as when normalization using the all‐genes approach was carried out (Fig. 3B).
Yeast stationary phase
To rule out the possibility that global changes are restricted to mammalian cell cultures, we also studied stationary phase in S. cerevisiae. This showed a 1.8‐fold median drop in mRNA levels when normalized using external controls (Fig. 4A). Significantly, there is a group of more than 1,000 genes that appear to be upregulated when normalized using the all‐genes method, although their mRNA levels had actually decreased. This was only revealed after normalization using external controls was carried out (Fig. 4C, left and middle columns). The apparent upregulation when the experiment was normalized using the all‐genes approach was seen because the amount of downregulation for these genes is less than that of most transcripts. When cellular RNA content was also taken into account, the changes were found to be more extreme (Fig. 4C, right column).
Changes dependent on expression levels
An important issue highlighted by the stationary phase experiment is that of an inconsistent distribution of changes across the range of expression levels: that is, along the x axis in Fig. 4A,B. This results in the apparently aberrant behaviour of the controls seen when normalization using the all‐genes approach was carried out (Fig. 4B). Highly expressed genes (on the right halves of the graphs in Fig. 4) show far more downregulation. One group of highly expressed mRNAs are the ribosomal protein genes (Holstege et al., 1998).The high degree of downregulation of these genes is consistent with the significant decrease in translation that occurs during stationary phase (Dickson & Brown, 1998). Unlike approaches using endogenous transcripts, normalization using external controls is not confounded by an uneven distribution of changes, and also contributes to a more accurate determination of mRNA changes in this respect.
Quantitative testing of reported changes
How accurate are the changes observed when normalization is carried out using external controls? Carrying out RT–PCR (reverse transcription followed by PCR) on a selected set of mRNAs leads to the question of what reference transcript should be chosen for the normalization of such experiments. The accuracy of the reported global mRNA changes was therefore tested by spiking the control mixtures in different ratios in identical pairs of RNA samples. The ratios chosen were 1:1, 1:2 and 1:10, thereby simulating 1‐, 2‐ and 10‐fold changes in the mRNA population relative to the controls. Normalization using external controls allowed the accurate determination of changes in mRNA levels relative to the controls. The gene spots drifted away from the control spots as the spiking ratio was increased to 1:10 (Fig. 5A–C). The median amount of change observed in the expression of the genes after normalization was almost identical to the spiking ratios, with values of 1.0‐, 1.9‐ and 9.3‐fold, respectively, for the average of each dye swap. This shows that the global changes reported here (up to twofold median changes; Figs 2,3,4) are likely to be correct.
This study focuses on the features used for single‐slide normalization and shows how the incorporation of external controls can markedly alter the interpretation of microarray experiments.
External controls have been used previously for normalization in studies examining the artificial inactivation of RNA polymerase II (Holstege et al., 1998; Wang et al., 2002). In these cases, global changes were expected. The experiments presented here show that global changes can also occur under more conventional experimental conditions. We have also shown how these controls can be applied to cope with local, intensity‐dependent systematic variation (Yang et al., 2002) by representation in sufficient numbers on each microarray subgrid, and by spiking over a range of levels. As well as being useful for normalization, controls such as these are useful for monitoring sample labelling, optimization of all microarray protocols, and as external controls for the reported ratios.
Recent papers have discussed the importance of experimental design and normalization algorithm choice (Kerr & Churchill, 2001; Quackenbush, 2001; Tseng et al., 2001; Kroll & Wolfl, 2002; Yang et al., 2002; Yang & Speed, 2002). An important improvement has been the adoption of normalization algorithms that take into account local, intensity‐dependent systematic variation (Yang et al., 2002). Regardless of the algorithm applied, most current analyses rely on assumptions of evenly distributed changes and/or on the absence of global shifts. As well as the wide use of the level of expression of all genes as an invariant feature, alternative normalization features that have been proposed include housekeeping genes, spotted microarray sample pools or spots containing genomic DNA (for overviews, see Kroll & Wolfl, 2002; Yang et al., 2002). Normalization using such features does not reveal global, unbalanced changes. The use of a common reference sample, made up of a collection of all the DNAs represented on arrays, has also been proposed (Dudley et al., 2002; Sterrenburg et al., 2002). This allows better comparisons between slides and also overcomes the problem of obtaining negligible signals for underrepresented mRNAs. However, such approaches do not address the possibility of global, unbalanced changes. It is likely that the use of a combination of external controls and a common pooled, spotted reference sample may be the best way of overcoming several problems simultaneously.
Monitoring global effects is not required, or feasible, for every microarray experiment. The advantage of the use of external controls is that it allows the possibility of detecting such changes in an experimental set‐up that is otherwise unchanged. The disadvantage is the requirement for the robust preparation of RNA samples and the reliable quantitation of yields. When RNA yields are too low to be monitored reliably, or when RNA preparations vary qualitatively, the use of external controls for normalization will not be reliable.
The importance of monitoring global and/or unbalanced mRNA changes is dictated by the goals of the experiment. For example, in disease classification studies, determination of the most extreme (relative) markers of a particular state without consideration of the actual changes does not require external control normalization. Neither is it important if the goal is to screen for only the most responsive genes in particular experimental conditions. However, as experimental biology becomes more comprehensive and quantitative due to the availability of whole‐genome sequences and to the increased focus on systems biology (Ideker et al., 2001), it will become more of a necessity to monitor with greater precision the actual, rather than the relative, changes in different cellular states. Methods that take into account the possibility of global changes will contribute towards such goals. Examples of studies that will benefit from the use of external controls include comprehensive studies of gene regulation and analyses of drug side‐effects, as well as microarray studies incorporating only limited sets of genes.
Accession numbers and protocols.
For MIAME (minimum information about a microarray experiment)‐compliant (Brazma et al., 2001) protocols, data sets in Microarray Gene Expression Markup Language (MAGE‐ML; Spellman et al., 2002) and normalization scripts, see our website (http://www.genomics.med.uu.nl/pub/jvp/ext_controls) or the public microarray database ArrayExpress (http://www.ebi.ac.uk/microarray/ArrayExpress/arrayexpress.html) (Table 1).
Constructs containing Bacillus subtilis genes (ycxA, yceG, ybdO, ybbR, ybaS, ybaF, ybaC, yacK and yabQ) cloned between the Xba I and Bam HI sites in pT7T3 (Amersham Pharmacia Biotech) were made, with an additional 30‐nucleotide poly(A) sequence § between the gene and the Xho I site. For making RNA, plasmids were digested with Xho I for use in in vitro transcription reactions using MEGAscript‐T7 (Ambion).
HUVECs were isolated as described in Jaffe et al. (1973). Cells were cultured in endothelial growth medium (EGM‐2) at 37 °C in the presence of 5% CO2. Before heat shock, the medium was removed, preheated EGM‐2 was added and cells were incubated at 42.5 °C in EGM‐2.
MCF7 cells were cultured in DMEM/Ham's F12 medium (1:1) containing 5% fetal calf serum, glutamine (300 mg ml−1), penicillin (100 inhibitory units ml−1) and streptomycin (100 mg ml−1). Cells at 70% confluence were grown for 30 h in phenol‐red‐free, serum‐free medium with 0.2% BSA, transferrin (10 mg ml−1) and 30 nM sodium selenite.
S. cerevisiae S288c (MATa; met15; ura3; his3Δ 1; leu2) (Research Genetics) was grown in YEP medium (containing yeast extract and peptone) supplemented with 2% glucose. Cultures for the spiking experiments (Fig. 5) were grown to mid‐log phase (OD600 = 0.5), and for the stationary phase experiment were grown to mid‐log phase or to stationary phase (OD600 = 10.0; 10‐day culture).
RNA isolation and labelling.
For mammalian cells, total RNA was prepared using Trizol (Gibco BRL) in accordance with the manufacturer's instructions. External controls were added in an appropriately diluted 5‐μl mixture to 10 μg of total RNA. Yeast total RNA was prepared using hot phenol. External controls were added, as an appropriately diluted 5‐μl mixture, to 500 μg of total RNA, and mRNA was isolated using Oligotex (Qiagen). Complementary DNA synthesis was carried out using 10 μg of mammalian total RNA, or 3 μg of yeast mRNA, in the presence of 2‐aminoallyl‐dUTP. Samples were purified using Microcon‐30 (Millipore) columns and were coupled to Cy3 and Cy5 fluorophores. Before hybridization, free dyes were removed using Chromaspin‐30 (Clontech) columns, and the efficiency of cDNA synthesis and dye incorporation was measured using a spectrophotometer (UV1240mini, Shimadzu).
From each sample, 300 ng cDNA (with a specific activity of 2–4% dye‐labelled nucleosides) was hybridized for 16–20 h at 42 °C. Slides were scanned in a Scanarray 4000 XL (Perkin Elmer Biosystems). Image analysis was carried out using Imagene 4.0 (Biodiscovery).
C6‐amino‐linked oligonucleotides (70 nucleotides in length), the Yeast Genome ArrayReady and the Human Genome ArrayReady Oligo set (version 1.1) were purchased from Qiagen, and were printed on Corning UltraGAPS slides with a MicroGrid II (Apogent Discoveries) using 48‐quill pins (Microspot2500; Apogent Discoveries) in 3 x SSC at 50% humidity and at 18 °C, and were processed by ultraviolet crosslinking (2,400 millijoules, 10 min) with a Stratalinker2400 (Stratagene).
Algorithms were based on Lowess print‐tip normalization (Yang et al., 2002), appliedin the statistical package R (Ihaka, 1996), using the existing packages SMA (http://www.stat.berkeley.edu/users/terry/zarray/Software/smacode.html) and com.braju.sma (http://www.maths.lth.se/help/R/com.braju.sma). Alterations were made for the import of Imagene 4.0 files, flagging of control spots, Lowess line calculation on subsets of spots (either controls or genes) and extrapolation to all spots in the subgrid. This has been incorporated in an R package called genomics.sma.
The first method, normalization using the expression levels of endogenous genes, uses gene spots to calculate the Lowess line for each subgrid and then applies these lines to all spots. The second method, normalization using the expression levels of all genes, and external controls, also shifts all spots linearly according to the median ratio of the control spots. The third method, normalization using external controls, uses the controls to calculate the Lowess line for each subgrid and then applies these lines to all spots for normalization.
The second and third methods gave almost identical results for the experiments shown in Figs 2, 3 and 5. Due to the non‐uniform distribution of changes in mRNA levels in the stationary phase experiment, only the third method (normalization using external controls) gave accurate results overall. This would be the method of choice for all experiments. However, there was a significant variation in the amount of the individual oligonucleotides supplied in the collections compared with the control oligonucleotides. This resulted in suboptimal extrapolation of the Lowess line to those gene spots that have extremely low intensity values, as a result of both low amounts of oligonucleotides and low or non‐existent levels of mRNA. The variation in the amount of oligonucleotides supplied has been rectified in later versions of the collections.
To take into account differences in RNA yield per cell, cells were counted, and a linear shift was applied to the normalized data using the ratio derived from the total RNA yield per cell from each pair of samples.
We thank A. Dijk, R. de Roos, and I. Hamelers for cell cultures, K. Duran for technical assistance, E. Holloway, H. Parkinson and U. Sarkans for assistance with MAGE‐ML generation and submission to ArrayExpress, C. Wilson and R.A. Young for external control constructs and P.C. van der Vliet, H.Th.M. Timmers, W. van Driel, R. Kerkhoven and C. Wijmenga for advice and discussions. This work was supported by the following grants from the Netherlands Organization for Scientific Research (NWO): 05050205, 90101238, 90101226, 016026009 and 90104219.
- Copyright © 2003 Nature Publishing Group