Evaluation of normalization methods for microarray data
© Park et al 2003
Received: 13 May 2003
Accepted: 02 September 2003
Published: 02 September 2003
Skip to main content
© Park et al 2003
Received: 13 May 2003
Accepted: 02 September 2003
Published: 02 September 2003
Microarray technology allows the monitoring of expression levels for thousands of genes simultaneously. This novel technique helps us to understand gene regulation as well as gene by gene interactions more systematically. In the microarray experiment, however, many undesirable systematic variations are observed. Even in replicated experiment, some variations are commonly observed. Normalization is the process of removing some sources of variation which affect the measured gene expression levels. Although a number of normalization methods have been proposed, it has been difficult to decide which methods perform best. Normalization plays an important role in the earlier stage of microarray data analysis. The subsequent analysis results are highly dependent on normalization.
In this paper, we use the variability among the replicated slides to compare performance of normalization methods. We also compare normalization methods with regard to bias and mean square error using simulated data.
Our results show that intensity-dependent normalization often performs better than global normalization methods, and that linear and nonlinear normalization methods perform similarly. These conclusions are based on analysis of 36 cDNA microarrays of 3,840 genes obtained in an experiment to search for changes in gene expression profiles during neuronal differentiation of cortical stem cells. Simulation studies confirm our findings.
Biological processes depend on complex interactions between many genes and gene products. To understand the role of a single gene or gene product in this network, many different types of information, such as genome-wide knowledge of gene expression, will be needed. Microarray technology is a useful tool to understand gene regulation and interactions [1–3]. For example, cDNA microarray technology allows the monitoring of expression levels for thousands of genes simultaneously. cDNA microarrays consist of thousands of individual DNA sequences printed in a high density array on a glass slide. After being reverse-transcribed into cDNA and labelled using red (Cy5) and green (Cy3) fluorescent dyes, two target mRNA samples are hybridized with the arrayed DNA sequences or probes. Then, the relative abundance of these spotted DNA sequences can be measured. After image analysis, for each gene the data consist of two fluorescence intensity measurements, (R, G), showing the expression level of the gene in the red and green labelled mRNA samples. The ratio of the fluorescence intensity for each spot represents the relative abundance of the corresponding DNA sequence. cDNA microarray technology has important applications in pharmaceutical and clinical research. By comparing gene expression in normal and tumor tissues, for example, microarrays may be used to identify tumor-related genes and targets for therapeutic drugs .
In microarray experiments, there are many sources of systematic variation. Normalization attempts to remove such variation which affects the measured gene expression levels. Yang et al.  and Yang et al.  summarized a number of normalization methods for dual labelled microarrays such as global normalization and locally weighted scatterplot smoothing (LOWESS ). Quackenbush  and Bilban et al.  provided good reviews on normalization methods. There have been some extensions for global and intensity-dependent normalizations. For example, Kepler et al.  considered a local regression to estimate a normalized intensities as well as intensity dependent error variance. Wang et al.  proposed a iterative normalization of cDNA microarray data for estimating a normalized coefficients and identifying control genes. Workman et al.  proposed a roust non-linear method for normalization using array signal distribution analysis and cubic splines. Chen et al.  proposed a subset normalization to adjust for location biases combined with global normalization for intensity biases. Edwards  considered a non-linear LOWESS normalization in one channel cDNA microarrays mainly for correcting spatial heterogeneity. The main idea of normalization for dual labelled arrays is to adjust for artifactual differences in intensity of the two labels. Such differences result from differences in affinity of the two labels for DNA, differences in amounts of sample and label used, differences in photomultiplier tube and laser voltage settings and differences in photon emission response to laser excitation. Although normalization alone cannot control all systematic variations, normalization plays an important role in the earlier stage of microarray data analysis because expression data can significantly vary from different normalization procedures. Subsequent analyses, such as differential expression testing would be more important such as clustering, and gene networks, though they are quite dependent on a choice of a normalization procedure [1, 3].
Several normalization methods have been proposed using statistical models (Kerr et al. ; Wolfinger et al. ). However, these approaches assume additive effects of random errors, which needs to be validated. Because they are less frequently used, we have not evaluated them here. Although several normalization methods have been proposed, no systematic comparison has been made for the performance of these methods. In this paper, we use the variability among the replicated slides to compare performance of several normalization methods.
We focus on comparing the normalization methods for cDNA microarrays. A detailed description on normalization methods considered in our study is given in the next section. Complex methods do not necessarily perform better than simpler methods. Complex methods may add noise to the normalized adjustment and may even add bias if the assumptions are incorrect. The fact that a non-linear method linearizes a graph of red intensity versus green intensity or a plot of log(R/G) versus log(RG) is not in itself evidence that the method is not just fitting noise. Consequently, proposed normalization methods require validation.
We first compare the methods using data from cDNA microarrays of 3,840 genes obtained in an experiment to search for changes in gene expression profiles during neuronal differentiation of cortical stem cells. Then, we perform simulation studies to compare these normalization methods systematically. Bolstad et al.  compare normalization methods for high density oligonucleotide array data using variance and bias.
The paper is organized as follows. Section 2 describes normalization methods. Section 3 describes the measures for variability to compare normalization methods. Section 4 shows the comparison results for cDNA microarrays obtained from a cortical stem cells experiment along with some simulation results. Finally, Section 5 summarizes the concluding remarks.
As pointed out by Yang et al. , the purpose of normalization is to remove systematic variation in a microarray experiment which affects the measured gene expression levels. They summarized a number of normalization methods: (i) within-slide normalization, (ii) paired-slide normalization for dye-swap experiments, and (iii) multiple slide normalization.
As a first step, one needs to decide which set of genes to use for normalization. Yang et al.  suggested three approaches: all genes on the array, constantly expressed genes, and controls. Recently, Tseng et al.  suggested using the rank invariant genes.
Global normalization (G) using the global median of log intensity ratios
Intensity dependent linear normalization (L)
Intensity dependent nonlinear normalization (N) using a LOWESS curve.
Under ideal experimental conditions, we expect M = 0 for the genes used for normalization. However, quite different patterns may be observed in real microarray experiments. Note that all three normalization methods consider regression models in the form of
For global normalization α0 is estimated by the median of M. For intensity dependent normalization (β0, β1) by least squares estimation, and c(A) by using the robust scatter plot smoother LOWESS. Thus, all three normalization methods are based on regression models of M in terms of A. Let (A) be the fitted value of c(A). The normalization process can be described in terms of (A). For gene j, the normalized ratio is
(A j ) = M j - (A j ). (1)
Abbreviation for Normalization Methods
Global median normalization
Intensity dependent linear normalization
Intensity dependent nonlinear normalization (LOWESS)
Print-tip scale normalization
Between-slide scale normalization
List of Normalization Methods
List of normalization methods including print-tip scale normalization
Global median normalization
Global median normalization on each print-tip
Global median normalization on each print-tip with scale normalization
Global median normalization and between-slide scale normalization
Global median normalization on each print-tip and between-slide scale normalization
Global median normalization on print-tip with scale normalization and between-slide scale normalization
Intensity dependent linear regression normalization
Intensity dependent linear regression normalization on each print-tip
Intensity dependent linear regression normalization on each print-tip with scale normalization
Intensity dependent linear regression normalization and between-slide scale normalization
Intensity dependent linear regression normalization on each print-tip and between-slide scale normalization
Intensity dependent linear regression normalization on each print-tip with scale normalization and between-slide scale normalization
Intensity dependent nonlinear regression normalization (LOWESS)
Intensity dependent nonlinear regression normalization (LOWESS) on each print-tip
Intensity dependent nonlinear regression normalization (LOWESS) on each print-tip with scale normalization
Intensity dependent nonlinear regression normalization (LOWESS) and between-slide scale normalization
Intensity dependent nonlinear regression normalization (LOWESS) on each print-tip and between-slide scale normalization
Intensity dependent nonlinear regression normalization (LOWESS) on each print-tip with scale normalization and between-slide scale normalization
In order to derive measures of variation, we now use y instead of M to represent the logarithm of the ratio of red and green background-corrected intensities. We introduce subscripts i,j,k and I for describing the cDNA microarray data introduced in the next section. Suppose there are I experimental groups denoted by i, J time points denoted by j, and K replications denoted by k.
Also suppose that there are N genes in each slide. For the simple case when there are I experimental groups but no time sequences, we can simply let J = 1. For the case when there are only time sequences, we can let I = 1. Thus, these subscripts represent general cases.
Let y ijkl be the logarithm of the red to green intensity ratio from group i(= 1, ..., I), time j(= 1, ... , J), replication k(= 1, ..., K), and gene l. Using the replicated observations, we want to derive the variability measures of y ijkl . It is expected that the better the normalization method, the smaller the variation among the replicated observations. Let σ l be the variability measure for the lth gene. We introduce two methods for estimating , l = 1, ... , N. Either density plots or box plots for the variability measures can be used for visually comparing different normalization methods.
Method 1. Pooled variance estimators
For gene l, a simple variance estimator for can be obtained by pooling variance estimators for each i and j. That is,
where , for i = 1, ..., I; j = 1, ..., J; l = 1, ..., N. Then, we compare the distributions of . The better the normalization method, the smaller the variance estimates.
Method 2. Variance estimator using analysis of variance models
Consider the following two-way analysis of variance(ANOVA) model with interactions for each gene.
y ijkl = μ l + α il + β jl + (αβ) ijl + ε ijkl , (3)
where i = 1, ..., I, j = 1, ..., J, k = 1, ..., K, and l = 1, ..., N. The gene effects μ l capture the overall mean intensity in fluorescent signals for genes across the arrays, groups, and time points. The α il terms account for gene specific group effects representing overall differences between two groups. The β jl account for time effects that capture differences in the overall concentration of mRNA in the samples from the different time points. The terms (αβ) ijl account for the interaction effects between group and time representing the signal contribution due to the combination of group and time.
From this ANOVA model, an unbiased estimate of can be obtained which is the error sum of squares divided by an appropriate degrees of freedom. Let be the variance estimate for the lth gene. Note that for Model (3) the values of are the same with those of (2), if there are no missing observations. However, this measure of variability is more flexible to use in the sense that it can be used with any ANOVA model which fits the intensity data well.
For this dataset, Park et al.  presented clustering analysis results after nonlinear LOWESS normalization. We applied three normalization methods to all 36 slides and compared the performance of normalization methods using the measures of variation among the replicated observations described in Section 3. This data set is a special case of the one described in Section 3 with two experimental groups, I = 2, six time points, J = 6, and three replications K = 3. For these microarray data, we derived the distribution of variance estimates by the methods given in Section 3. Since the two methods provided quite similar results, we only present the results from the ANOVA model.
As shown in Figure 3a, the three normalization methods G, L, N reduce the variability greatly. The two intensity dependent normalization methods denoted by L and N perform equally well and somewhat better than global normalization.
Figure 3b shows the results of global normalization with different scale normalization approaches. Great differences are not found among the six approaches. Hence global normalization does not appear to be improved by normalizing separately over the grids of spots defined by the print-tips or by scale adjustments.
Figure 3c focuses on intensity dependent linear normalization with different scale normalization approaches. Figure 3d is for intensity dependent non-linear LOWESS normalization with different scale normalization approaches. Neither linear nor non-linear intensity dependent normalization appear to benefit greatly by scale adjustments or by performing normalization by print-tip grids for this set of data.
Based on these figures, we conclude that intensity dependent normalization methods perform better than global normalization methods. Small differences are observed between the linear and nonlinear normalization methods. In addition, small effects are observed for scale normalization methods.
Figure 4a shows the plot of (logG, logR) for Type I which shows a clear linear pattern. Most observations are scattered around the y = x line showing a high correlation between the two intensity values.
Figure 4b shows the second type. It has a specific non-linear pattern with many spots distant from the y = x line.
Figure 4c shows the third type which has the same linear pattern as Type I. However, the lower intensity values have higher variabilities than those of higher intensity values.
Figure 4d shows the fourth type. It is a combination of Type II and Type III. That is, it has a non-linear pattern with higher variability at the lower intensities. In order to generate (logG I , logR I ), we assume that they are observed from the bivariate normal distribution with mean vector β and covariance matrix Σ, where
For simplicity, let the random vector (Y G , Y R ) T represent the log-transformed intensity values, say (logG I , logR I ). To generate bivariate normal random variables, it is convenient to generate Y G first and then generate Y R from the conditional distribution of Y R | Y G = Y. Note that the conditional distribution is given by
where μ R|G = μ R - ρσ R /σ G (Y - μ G ) and . For Type II, we use the same values of Y G generated from Type I. For Y R, however, we use the transformed mean and variance to add a pattern like that shown in Figure 4b. That is, and . We first tried generating Type II data by directly transforming the data generated from Type I but could not get a graph with the desired pattern. The functions f 1 and f 2 are found by trial and error. For Type III, we use a similar approach. That is, we use the same values of Y G generated from Type I. For Y R , we use the transformed variance to allow larger variability for the lower intensity observations. Thus, and . For simplicity, we use the same function f 2 used for Type II. Similarly, for Type IV we use the same values of Y G generated from Type I. For Y R , we use the transformed mean and variance to add a specific pattern. That is, and . The functions f 3 and f 4 are also found by trial and error. For each type, we generate ten replicated samples. For these replicated samples, we apply the normalization methods and draw the distributions for the variation measures. For simplicity, we do not consider the print-tip variation in this simulation. We compare four graphs: original data(O), globally normalized data(G), linear normalized data(L), and non-linear normalized data(N) using LOWESS.
We summarize our findings from the simulation studies. First, the normalization methods performed similarly when the original data has a linear pattern such as Type I, and when there are high variabilities for observations with lower intensities such as Types III. Second, when there is a specific pattern such as Type II, all normalization methods tend to reduce MSEs. Third, in this case, MSEs are dominated by biases. Fourth, only small differences are observed between the intensity-dependent linear and non-linear normalization methods.
In this article, we compare normalization methods commonly used to analyze microarray data. The comparison is based on the variability measures derived from the replicated microarray samples. These variability measures can be easily derived from any replicated microarray experiment. As pointed out by many researchers, many undesirable systematic variations are observed in the microarray experiment. Normalization becomes a standard process for removing some of the variation which affects the measured gene expression levels. Although a number of normalization methods have been proposed, it has been difficult to decide which method performs better than the others. Thus, the evaluation of normalization methods in microarray data analysis is indeed an important issue. In this article, we show that the intensity dependent normalization method performs better than the simpler global normalization methods in many cases. We have not been sure about whether apparent nonlinearity of an M-A scatter plot or a scatter plot of red vs green is sufficient basis for feeling confident that a non-linear normalization is useful. Although we have studied only a limited number of data sets, our findings can provide some guidance on the selection of normalization methods. There are clearly cases where intensity dependent normalization performs substantially better than global normalization methods. In most of the cases that we considered, the non-linear intensity dependent procedures did not perform substantially better than a linear intensity dependent method, although there may be datasets where this is not the case. For the cases we considered, we did not see large benefits to separate normalization by print-tip grids or for scale normalization. None of the normalization methods effectively addressing the dependence of the variance of measurements on intensity level. Recently, some other non-linear normalization methods have been employed such as B-splines and Gaussian-kernel fitting [12, 21]. These non-linear normalization methods seem to be comparable to the LOWESS normalization. We think the basic idea of non-linear normalization is the same whether they use LOWESS curve, splines, or kernels. They are quite effective in controlling specific non-linear variations. Some methods may have more computational efficiency. We believe that more sophisticated non-linear normalization methods are under development. We recommend that experimentalists examine their data carefully and consider applying intensity-dependent normalization methods routinely. Software for intensity dependent normalization is available on desktop packages such as BRB-ArrayTools http://linus.nci.nih.gov/BRB-ArrayTools and via the internet in SNOMAD http://pevsnerlab.kennedykrieger.org/snomad.htm tools which provide any DNA array researcher with the tools to apply a variety of non-linear intensity-dependent normalization methods to their data .
The authors wish to thank two anonymous referees whose comments were extremely helpful. The work was supported by the IMT2000 contribution from Korea Ministry of Health and Welfare, and Ministry of Information and Communication (01-PJ11-PG9-01BT-00A-0045).
This article is published under license to BioMed Central Ltd. This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL.