Access the full text.
Sign up today, get DeepDyve free for 14 days.
E. Schadt, C. Li, C. Su, W. Wong (2001)
Analyzing high‐density oligonucleotide gene expression array dataJournal of Cellular Biochemistry, 80
(2001)
Probe level quantile normalization of high density oligonucleotide array data. Unpublished Manuscript
(2000)
Large-scale genomic analysis using Affymetrix GeneChip R
E. Schadt, C. Li, B. Ellis, W. Wong (2001)
Feature extraction and normalization algorithms for high‐density oligonucleotide gene expression array dataJournal of Cellular Biochemistry, 84
R. Lipshutz, S. Fodor, T. Gingeras, D. Lockhart (1999)
High density synthetic oligonucleotide arraysNature Genetics, 21 Suppl 1
R. Irizarry, B. Hobbs, F. Collin, Y. Beazer-Barclay, Kristen Antonellis, U. Scherf, T. Speed (2003)
Exploration, normalization, and summaries of high density oligonucleotide array probe level data.Biostatistics, 4 2
Cheng Li, Wing Wong (2001)
Model-based analysis of oligonucleotide arrays: model validation, design issues and standard error applicationGenome Biology, 2
W. Cleveland, S. Devlin (1988)
Locally Weighted Regression: An Approach to Regression Analysis by Local FittingJournal of the American Statistical Association, 83
A. Hartemink, D. Gifford, T. Jaakkola, R. Young (2001)
Maximum-likelihood estimation of optimal scaling factors for expression array normalization, 4266
W. Venables, B. Ripley (1996)
Modern Applied Statistics with S-Plus.Biometrics, 52
S. Dudoit, Y. Yang, M. Callow, T. Speed (2002)
STATISTICAL METHODS FOR IDENTIFYING DIFFERENTIALLY EXPRESSED GENES IN REPLICATED cDNA MICROARRAY EXPERIMENTSStatistica Sinica, 12
W. Venables, B. Ripley (1997)
Modern Applied Statistics with S-Plus Second edition
(2001)
Statistical algorithms reference guide
I. Sidorov, D. Gee, D. Dimitrov, Douglas Hosack, J. Yang, R. Lempicki, M. Cam (2002)
Oligonucleotide Microarray Data Distribution And NormalizationInf. Sci., 146
(2002)
GeneLogic
Ross Ihaka (1996)
Gentleman R: R: A language for data analysis and graphics
M. Schena (2000)
Microarray Biochip Technology
Vol. 19 no. 2 2003 BIOINFORMATICS Pages 185–193 A comparison of normalization methods for high density oligonucleotide array data based on variance and bias 1,∗ 2 3 4, 5 B. M. Bolstad ,R.A. Irizarry ,M. Astrand and T. P. Speed Group in Biostatistics, University of California, Berkeley, CA 94720, USA, Department of Biostatistics, John Hopkins University, Baltimore, MD, USA, 3 4 AstraZeneca R & D Molndal, ¨ Sweden, Department of Statistics, University of California, Berkeley, CA 94720, USA and Division of Genetics and Bioinformatics, WEHI, Melbourne, Australia Received on June 13, 2002; revised on September 11, 2002; accepted on September 17, 2002 ABSTRACT oligonucleotides of 25 base pairs in length are used to Motivation: When running experiments that involve multi- probe genes. There are two types of probes: reference ple high density oligonucleotide arrays, it is important to re- probes that match a target sequence exactly, called the move sources of variation between arrays of non-biological perfect match (PM), and partner probes which differ from origin. Normalization is a process for reducing this varia- the reference probes only by a single base in the center of tion. It is common to see non-linear relations between ar- the sequence. These are called the mismatch (MM) probes. rays and the standard normalization provided by Affymetrix Typically 16–20 of these probe pairs, each interrogating a does not perform well in these situations. different part of the sequence for a gene, make up what Results: We present three methods of performing nor- is known as a probeset. Some more recent arrays, such malization at the probe intensity level. These methods are as the HG-U133 arrays, use as few as 11 probes in a called complete data methods because they make use of probeset. The intensity information from the values of data from all arrays in an experiment to form the normaliz- each of the probes in a probeset are combined together ing relation. These algorithms are compared to two meth- to get an expression measure, for example, Average ods that make use of a baseline array: a one number scal- Difference (AvgDiff), the Model Based Expression Index ing based algorithm and a method that uses a non-linear (MBEI) of Li and Wong (2001), the MAS 5.0 Statistical normalizing relation by comparing the variability and bias algorithm from Affymetrix (2001), and the Robust Multi- of an expression measure. Two publicly available datasets chip Average proposed in Irizarry et al. (2003). are used to carry out the comparisons. The simplest and The need for normalization arises naturally when quickest complete data method is found to perform favor- dealing with experiments involving multiple arrays. There ably. are two broad characterizations that could be used for the Availabilty: Software implementing all three of the com- type of variation one might expect to see when comparing plete data normalization methods is available as part of arrays: interesting variation and obscuring variation. the R package Affy, which is a part of the Bioconductor We would classify biological differences, for example project http://www.bioconductor.org. large differences in the expression level of particular Contact: bolstad@stat.berkeley.edu genes between a diseased and a normal tissue source, Supplementary information: Additional figures may be as interesting variation. However, observed expression found at http://www.stat.berkeley.edu/∼bolstad/normalize/ levels also include variation that is introduced during index.html the process of carrying out the experiment, which could be classified as obscuring variation. Examples of this INTRODUCTION obscuring variation arise due to differences in sample The high density oligonucleotide microarray technology, preparation (for instance labeling differences), production of the arrays and the processing of the arrays (for instance as provided by the Affymetrix GeneChip ,is being used scanner differences). The purpose of normalization is in many areas of biomedical research. As described in to deal with this obscuring variation. A more complete Lipshutz et al. (1999) and Warrington et al. (2000), discussion on the sources of this variation can be found in To whom correspondence should be addressed. Hartemink et al. (2001). Bioinformatics 19(2) c Oxford University Press 2003; all rights reserved. 185 B.M.Bolstad et al. Affymetrix has approached the normalization problem M = log x /x and A = log x x .A k ki kj k ki kj 2 2 by proposing that intensities should be scaled so that each normalization curve is fitted to this M versus A plot array has the same average value. The Affymetrix nor- using loess. Loess is a method of local regression (see malization is performed on expression summary values. Cleveland and Devlin 1988 for details). The fits based on This approach does not deal particularly well with cases the normalization curve are M and thus the normalization where there are non-linear relationships between arrays. adjustment is M = M − M . Adjusted probe intensites k k Approaches using non-linear smooth curves have been M M K k A + A − k K 2 2 are given by x = 2 and x = 2 . proposed in Schadt et al. (2001, 2002) and Li and Wong ki kj (2001). Another approach is to transform the data so that The preferred method is to compute the normalization the distribution of probe intensities is the same across a curves using rank invariant sets of probes. This paper uses set of arrays. Sidorov et al. (2002) propose parametric invariants sets since it increases the implementation speed. and non-parametric methods to achieve this. All these To deal with more than two arrays, the method is approaches depend on the choice of a baseline array. extended to look at all distinct pairwise combinations. The We propose three different methods of normalizing normalizations are carried out in a pairwise manner as probe intensity level oligonucleotide data, none of which above. We record an adjustment for each of the two arrays is dependent on the choice of a baseline array. Normaliza- in each pair. So after looking at all pairs of arrays for any tion is carried out at probe level for all the probes on an array k where 1 ≤ k ≤ n,wehave adjustments for chip k array. Typically we do not treat PM and MM separately, relative to arrays 1,..., k − 1, k + 1,..., n.We weight the but instead consider them all as intensities that need to be adjustments equally and apply to the set of arrays. We have normalized. The normalization methods do not account found that after only 1 or 2 complete iterations through all for saturation. We consider this a separate problem to be pairwise combinations the changes to be applied become dealt with in a different manner. small. However, because this method works in a pairwise In this paper, we compare the performance of our three manner, it is somewhat time consuming. proposed complete data methods. These methods are then Contrast based method The contrast based method is compared with two methods making use of a baseline array. The first method, which we shall refer to as the another extension of the M versus A method. Full details scaling method, mimics the Affymetrix approach. The can be found in Astrand (2001). The normalization is second method, which we call the non-linear method, carried out by placing the data on a log-scale and mimics the approaches of Schadt et al. Our assessment transforming the basis. In the transformed basis, a series of the normalization procedures is based on empirical of n − 1 normalizing curves are fit in a similar manner to results demonstrating ability to reduce variance without the M versus A approach of the cyclic loess method. The increasing bias. data is then adjusted by using a smooth transformation which adjusts the normalization curve so that it lies NORMALIZATION ALGORITHMS along the horizontal. Data in the normalized state is obtained by transforming back to the original basis and Complete data methods exponentiating. The contrast based method is faster than The complete data methods combine information from all the cyclic method. However, the computation of the loess arrays to form the normalization relation. The first two smoothers is still somewhat time consuming. methods, cyclic loess and contrast, are extensions of ac- cepted normalization methods that have been used suc- Quantile normalization The goal of the quantile method cessfully with cDNA microarray data. The third method, is to make the distribution of probe intensities for each based on quantiles, is both quicker and simpler than those array in a set of arrays the same. The method is motivated methods. by the idea that a quantile–quantile plot shows that the distribution of two data vectors is the same if the plot is Cyclic loess This approach is based upon the idea of a straight diagonal line and not the same if it is other than the M versus A plot, where M is the difference in a diagonal line. This concept is extended to n dimensions log expression values and A is the average of the log so that if all n data vectors have the same distribution, expression values, presented in Dudoit et al. (2002). then plotting the quantiles in n dimensions gives a straight However, rather than being applied to two color channels 1 1 √ √ line along the line given by the unit vector ,..., . on the same array, as is done in the cDNA case, it is n n applied to probe intensities from two arrays at a time. An This suggests we could make a set of data have the same M versus A plot for normalized data should show a point distribution if we project the points of our n dimensional cloud scattered about the M = 0 axis. quantile plot onto the diagonal. Forany two arrays i , j with probe intensities x and x Let q = (q ,..., q ) for k = 1,..., p be the vector ki kj k1 kn where k = 1,..., p represents the probe, we calculate of the kth quantiles for all n arrays q = (q ,..., q ) k1 kn 186 Normalization of Oligonucleotide Arrays intensities). Then the intensities for the normalized array 1 1 √ √ and d = ,..., be the unit diagonal. To transform n n would be from the quantiles so that they all lie along the diagonal, x = β x i i consider the projection of q onto d One can also easily implement the scaling algorithm by n n 1 1 using probes from a subset of probesets chosen by using proj q = q ,..., q kj kj some stability criteria. The HG-U133 arrays provide a set n n j =1 j =1 of probesets that have been selected for stability across tissue types, and these could be used for establishing a This implies that we can give each array the same normalization. distribution by taking the mean quantile and substituting it as the value of the data item in the original dataset. This Non-linear method The scaling method is equivalent to motivates the following algorithm for normalizing a set of fitting a linear relationship with zero intercept between the data vectors by giving them the same distribution: baseline array and each of the arrays to be normalized. This normalizing relation is then used to map from each 1. given n arrays of length p, form X of dimension array to the baseline array. This idea can be extended to use p × n where each array is a column; a non-linear relationship to map between each array and 2. sort each column of X to give X ; the baseline array. Such an approach is detailed in Schadt sort et al. (2002). This method is used in Li and Wong (2001) 3. take the means across rows of X and assign this sort and implemented in the dChip software http://www.dchip. mean to each element in the row to get X ; sort org. The general approach of these papers is to select a 4. get X by rearranging each column of normalized set of approximately rank invariant probes (between the X to have the same ordering as original X baseline and arrays to be normalized) and fit a non-linear sort relation, like smoothing splines as in Schadt et al. (2002), The quantile normalization method is a specific case of or a piecewise running median line as in Li and Wong −1 the transformation x = F (G (x )), where we estimate (2001). G by the empirical distribution of each array and F The non-linear method used in this paper is as follows. using the empirical distribution of the averaged sample First we select a set of probes for which the ranks are quantiles. Extensions of the method could be implemented invariant across all the arrays to be normalized. Then we −1 where F and G are more smoothly estimated. fit loess smoothers to relate the baseline to each of the One possible problem with this method is that it forces arrays to be normalized. These loess normalization curves the values of quantiles to be equal. This would be are then used to map probe intensities from the arrays to most problematic in the tails where it is possible that a be normalized to the baseline. This approach is intended probe could have the same value across all the arrays. to mimic the approach used in dChip. We expect loess However, in practice, since probeset expression measures smoothers to perform in the same manner as splines or are typically computed using the value of multiple probes, a running median line. we have not found this to be a problem. Suppose that f (x ) is the loess smoother mapping from array i to the baseline. Then, in the same notation as above, Methods using a baseline array the normalized array probe intensities are Scaling methods The standard Affymetrix normalization is a scaling method that is carried out on probeset x = f (x ) i i expression measures. To allow consistent comparison with our other methods, we have carried out a similar Note that as with the scaling method, the baseline is the normalization at the probe level. Our version of this array having the median of the median probe intensities. method is to choose a baseline array, in particular, the array having the median of the median intensities. All DATA arrays are then normalized to this ‘baseline’ via the We make use of data from two sets of experiments: following method. If x are the intensities of the base Adilution/mixture experiment and an experiment using baseline array and x is any array, then let spike-ins. We use these datasets because they allow us to x˜ assess bias and variance. The dilution/mixture and spike- base β = in datasets are available directly from GeneLogic (2002) x˜ and have been made available for public comparison of where x˜ is the trimmed mean intensity (in our analysis analysis methods. This data has been previously described we have excluded the highest and lowest 2% of probe in Irizarry et al. (2003). 187 B.M.Bolstad et al. Density of PM probe intensities for SpikeIn chips After Quantile Normalization 46 8 10 12 14 log(PM) Fig. 1. A plot of the densities for PM for each of the 27 Fig. 2. 10 pairwise M versus A plots using liver (at concentration spike-in datasets, with distribution after quantile normalization 10) dilution series data for unadjusted data. superimposed. Dilution/mixture data The dilution/mixture data series consists of 75 HG-U95A (version 2) arrays, where two sources of RNA, liver (source A), and a central nervous system cell line (source B) are investigated. There are 30 arrays for each source, broken into 6 groups at 5 dilution levels. The remaining 15 arrays, broken into 3 groups of 5 chips, involve mixtures of the two tissue lines in the following proportions: 75 : 25, 50 : 50, and 25 : 75. Spike-in data The spike-in data series consists of 98 HG-U95A (version 1) arrays where 11 different cRNA fragments have been spiked in at various concentrations. There is a dilution series consisting of 27 arrays which we will examine in Fig. 3. 10 pairwise M versus A plots using liver (at concentration this paper. The remaining arrays are two sets of latin 10) dilution series data after quantile normalization. square experiments, where in most cases three replicate arrays have been used for each combination of spike- in concentrations. We make use of 6 arrays (two sets of between arrays. The same 10 pairwise comparisons can triplicates) from one of the latin squares. be seen after quantile normalization in Figure 3. The point clouds are all centered around M = 0. Plots produced us- RESULTS ing the contrast and cyclic loess normalizations are similar. Probe level analysis Expression measures Figure 1 plots the densities for the log(PM ) for each of Comparing normalization methods at the probeset level the 27 arrays from the spike-in dataset, along with the requires that one must decide on an expression measure. distribution obtained after quantile normalization. Although in this paper we focus only on one expression An M versus A plot allows us to discern intensity de- pendent differences between two arrays. Figure 2 shows measure, the results obtained are similar when using other M versus A plots for unadjusted PM for all 10 possible measures. pairs of 5 arrays in the liver 10 group before normaliza- The expression summary used in this paper is a robust tion. Clear differences between the arrays can be seen by combination of background adjusted PM intensities and looking at the loess lines. The point clouds are not cen- is outlined in Irizarry et al. (2003). We call this method tered around M = 0 and we see non-linear relationships the Robust Multichip Average (RMA). RMA estimates Density 0.0 0.2 0.4 0.6 0.8 1.0 Normalization of Oligonucleotide Arrays Quantile:unnormalized Loess:Quantile are based upon a robust average of log B PM , where ( ( )) B (PM ) are background corrected PM intensities. The expression measure may be used on either the natural or log scales. Irizarry et al. (2003) contains a more complete discus- sion of the RMA measure, and further papers exploring its 24 6 8 10 12 14 2 468 10 12 14 properties are under preparation. log mean log mean Contrast:Quantile Contrast:Loess Probeset measure comparisons Variance comparisons In the context of the dilution study, consider the five arrays from a single RNA source within a particular dilution level. We calculate expression measures for every probeset on each array and then compute the variance and mean of the probeset expression 24 6 8 10 12 14 2 4 6 8 10 12 14 summary across the five arrays. This is repeated for each log mean log mean group of 5 arrays for the entire dilution/mixture study. We do this after normalization by each of our three complete Fig. 4. log variance ratio versus average log mean for liver 2 2 data methods. dilution data at concentration 10. Plotting the log of the ratio of variances versus the av- erage of the log of the mean (expression measure across arrays) allows us to see differences in the between array Quantile:Scaling variations and intensity dependent trends when comparing normalization methods. In this case, the expression mea- sures have all been calculated on the natural scale. Fig- ure 4 shows such plots for the liver at the dilution level 10. Specifically, the four plots compare the variance ratios 246 8 10 12 14 for quantile : unnormalized, loess : quantile, contrast : log mean quantile and contrast : loess. The horizontal line indicates the x -axis. The other line is a loess smoother. Where the Quantile:Nonlinear loess smoother is below the x -axis, the first method in the ratio has the smaller ratio and vice versa when the loess smoother is above the line. All three methods reduce the variance at all intensity levels in comparison to data that has not been normalized. The three normalization methods 468 10 12 14 perform in a relatively comparable manner, but the quan- log mean tile method performs slightly better for this dataset, as can be seen in the loess : quantile and contrast : quantile plots. Fig. 5. log variance ratio versus average log mean using the spike- Similar results are seen in comparable plots (not shown) 2 2 in data. Comparing the baseline methods with the quantile method. for the other dilution/mixture groups. We repeat this analysis with the 27 spike-in arrays, but this time we include the two baseline methods in our comparison. The complete data methods generally leave (or lower) due to the shifting and not because of the the mean level of a particular probeset at a level similar normalization. To minimize this problem and make a fairer to that achieved when using unnormalized data. However, comparison, we work with the expression measure on the when one of the two baseline methods is used, the mean log scale when comparing the baseline methods to the of a particular probeset is more reminiscent of the value of that probeset in the baseline array. In the natural scale, it is complete data methods. Figure 5 compares the baseline easy to see a mean-variance relationship, where a higher methods to the quantile methods. We see that the quantile mean implies high variability. Thus, when a comparison is method reduces the between array variances more than the made between the baseline methods and the complete data scaling method. The non-linear normalization performs methods, we find that if a baseline array which shifts the a great deal closer to the quantile method. Similar plots intensities higher (or lower) than the level of those of the (not shown) comparing the complete data methods with unnormalized means is selected, then the corresponding each other for the spike-in data demonstrate that quantile variance of the probeset measures across arrays is higher normalization has a slight edge over all the other methods. log variance ratio log variance ratio log variance ratio log variance ratio –4 –2 0 2 4 –10 –8 –6 –4–2 0 2 4 –4 –2 0 2 4 –4 –2 0 2 4 log variance ratio log variance ratio –4 –2 0 2 4 –4 –2 0 2 4 B.M.Bolstad et al. Spike in data (based on 352 pairwise comparisions) spike-in probesets, we fit the following linear model Quantile log E = β + β log c + 0 1 2 2 Contrast Loess Scaling where E is the value of the expression measure and c are nonlinear the concentrations. Note that the array with spike-in con- centration 0 is excluded from the model fit, although it is used in the normalization. The ideal results would be to have slopes that are near 1. Table 1 shows the slope esti- mates for each of the spike-in probesets after normaliza- tion by each of the three complete data methods, the two methods using a baseline and when no normalization has taken place. For the three complete data methods for 10 out of the 11 spike-in probesets, the quantile method gives 24 6 8 10 12 a slope closer to 1 and the non-linear method has slopes lower than the complete data methods. However, both the scaling and not normalizing have slopes closer to 1. For Fig. 6. Comparing the ability of methods to reduce pairwise the non-spike-in probesets on these arrays we should see differences between arrays by using average absolute distance from no linear relation if we fit the same model, since there loess smoother to x -axis in pairwise M versus A plots using spike-in should be no relation between the spike-in concentrations dataset. Smaller distances are favorable. and the probeset measures. Fitting the linear model above, we find that there is a median slope of 0.042 for these probesets using the unnormalized data. For the quantile A similar plot (not shown) comparing the two baseline method the value of the median slope is −0.005. All the methods shows as expected that the non-linear method other normalization methods have median slope near 0. reduces variance when compared to the scaling method. This is about the same difference in slopes as we observed for the spike-ins when comparing the unnormalized data Pairwise comparison The ability to minimize dif- and the best of the normalization methods. In other words, ferences in pairwise comparisons between arrays is a there is a systematic trend due to the manner in which the desirable feature of a normalization procedure. An M arrays were produced that has resulted in the intensities of versus A plot comparing expression measures on two all the probesets being related to the concentration of the arrays should be centered around M = 0if there is no spike-ins. We should adjust the spike-in slopes by these clear trend towards one of the arrays. Looking at the amounts. For example, we could adjust the slope of BioB- absolute distance between a loess for the M versus A 5 for the quantile method to 0.845 + 0.005 = 0.850 and plot and the x -axis allows us to assess the difference in the unnormalized slope to 0.893 − 0.042 = 0.851. array to array comparisons. We can compare methods by The average R for the spike-in probesets, excluding looking at this distance across a range of intensities and CreX-3, are 0.87 for the quantile method, and 0.855, averaging the distance across all pairwise comparisions. 0.849, 0.857 and 0.859 for the contrast, loess, non-linear Figure 6 shows such a plot for the spike-in data, we and scaling methods, respectively. It was 0.831 for the see that the scaling method performs quite poorly when unnormalized data. The median standard error for the compared to the three complete data methods. The non- slopes was 0.063 for the quantile method. For the other linear method performs at a similar level to the complete methods these standard errors were 0.065 (contrast), 0.068 data methods. For this dataset the quantile method is (cyclic loess), 0.063 (non-linear), 0.065 (scaling) and slightly better. An important property of the quantile 0.076 (unnormalized). Thus of all the algorithms, the method is that these differences remain relatively constant quantile method has high slopes, a better fitting model and across intensities. more precise slope estimates. Bias comparisons One way to look at bias is in the The slopes may not reach 1 for several reasons. It is context of the spike-in dilution series. We use data for the possible that there is a ‘pipette’ effect. In other words we 27 arrays from the spike-in experiment with 11 control can not be completely sure of the concentrations. It is more fragments spiked in at 13 different concentrations (0.00, likely that we observe concentration plus an error which 0.50, 0.75, 1.00, 1.50, 2.00, 3.00, 5.00, 12.50, 25.00, leads to a downward bias in the slope estimates. Other 50.00, 75.00, 100.00, 150.00 pM). We normalize each possible reasons include the saturation of signal at the high of the 27 arrays as a group using each of the quantile, end (this is not a concern with this data) and having a contrast, cyclic loess and scaling normalizations. To the higher background effect at the lower end. avg abs distance from lowess to x axis 0.0 0.1 0.2 0.3 0.4 Normalization of Oligonucleotide Arrays Table 1. Regression slope estimates for spike-in probesets. A slope closer to one is better Name Quantile Contrast Loess Non-linear Scaling None AFFX-BioB-5 at 0.845 0.837 0.834 0.803 0.850 0.893 AFFX-DapX-M at 0.778 0.771 0.770 0.746 0.783 0.826 AFFX-DapX-5 at 0.754 0.747 0.728 0.731 0.764 0.807 AFFX-CreX-5 at 0.903 0.897 0.889 0.875 0.912 0.955 AFFX-BioB-3 at 0.836 0.834 0.825 0.807 0.848 0.890 AFFX-BioB-M at 0.789 0.782 0.781 0.762 0.797 0.838 AFFX-BioDn-3 at 0.547 0.543 0.550 0.514 0.553 0.595 AFFX-BioC-5 at 0.801 0.794 0.793 0.763 0.808 0.851 AFFX-BioC-3 at 0.796 0.790 0.785 0.769 0.805 0.847 AFFX-DapX-3 at 0.812 0.804 0.793 0.776 0.815 0.859 AFFX-CreX-3 at −0.007 −0.006 0.002 −0.007 0.005 0.046 Non-spike-in (median) −0.005 −0.005 −0.005 −0.007 −0.001 0.042 The problems with choosing a baseline The non-linear Mean of probeset expression based on 6 chips (and scaling) method requires the choice of a baseline. In this paper we have chosen the array having the median median, but other options are certainly possible. To address these concerns we examine a set of six arrays chosen from the spike-in datasets. In particular, we choose two sets of triplicates, where the fold change of each of the spike-in probesets between the two triplicates is large. The two triplicates are chosen so that about half of the probesets are high in one triplicate and low in the other and vice versa. We normalize this dataset using both the quantile and non-linear methods. However, for the non-linear method we experiment with the use of each of the six arrays as the baseline array. We also try using two synthetic baseline None Quantile means median 1 2345 6 arrays: one constructed by taking probewise means Method and one taking probewise medians. Figure 7 shows the distribution of the mean of the probeset measure across Fig. 7. Distribution of average (over 6 chips) of a probeset arrays. We see that the quantile normalization produces expression measure using different baseline normalizations. a set of means that is very similar in distribution to the means of the unnormalized data. The means from the non-linear normalizations using each of the six different the baselines. However, two baseline arrays perform quite baseline arrays are quite different from each other and the poorly. As noted before, this is a reflection of the baseline unnormalized data. This is somewhat of a drawback to methods shifting the intensities higher or lower depending the baseline methods. It seems more representative of the on the baseline. The mean based synthetic baseline does complete data to consider all arrays in the normalization not reduce the variabilty of the probeset measure to the rather than to use only a single baseline and give the normalized data characteristics closer to those of one par- same degree as the quantile method. ticular array. Only the mean based synthetic baseline array Looking at the 11 spike-in probesets, we calculate bias comes close to the unnormalized and quantile methods. by taking the difference between the log of the ratio Table 2 summarizes some results from this analysis. between spike-in concentrations and the log of the ratio We see that all the methods reduce the variability of of intensities in the two groups. One spike-in, Crex-3, did the probeset measure between arrays compared to that not seem to perform quite as well as the other spike-ins of the unnormalized data. In each case, around 95% of and was excluded from the analysis. Looking at the total the probesets have reduced variance. When compared to absolute bias across the 10 spike-in probesets, we see that the quantile method, it comes out more even with a little the non-linear method has lower total bias (compared to over 50% of probesets having reduced variance for four of the quantile method) for four of the methods, but two 2468 10 12 14 B.M.Bolstad et al. Table 2. Comparing variance and bias with the non-linear normalization when using different baselines Method % with lower var % lower var Abs # abs # abs reduced cf. U reduced cf. Q Bias Bias cf U Bias cf Q Probewise mean 83 40 9.2 5 5 Probewise median 96 58 7.9 6 6 Non-linear 1 96 53 7.5 7 5 Non-linear 2 93 31 11.8 2 4 Non-linear 3 94 37 10.5 4 4 Non-linear 4 95 47 7.4 6 5 Non-linear 5 96 55 7.4 7 5 Non-linear 6 96 55 7.5 7 5 Quantile (Q) 95 NA 8.5 6 NA Unnormalized (U) NA NA 9.7 NA NA M vs A plot for spikein triplicates are even bigger than for unnormalized data. Again, this is related to the mean–variance relationship. The four baselines shifted slightly lower in the intensity scale give the most precise estimates. Using this logic, one could argue that choosing the array with the smallest spread and centered at the lowest level would be the best, but this does not seem to treat the data on all arrays fairly. Compared to the unnormalized data, 6 of the spike-in probesets from the quantile normalized data have a smaller bias. For the non- linear normalization using array 1 as the baseline (this is the array chosen using our heuristic), 7 had smaller bias. However, looking at the other baselines, anywhere from 2to7 probesets had lower bias. When compared to the quantile method, the results are more even, with about 2468 10 12 an equal number of the spike-in probesets having a lower bias when using the non-linear method as when using the quantile normalization. An M versus A plot between the Fig. 8. M versus A plot for spike-in triplicate data normalized using two groups shows all the spike-in points clearly outside the quantile normalization. Spike-ins are clearly identified. point cloud, no matter which normalization is used. This plot for quantile normalized data is shown in Figure 8. performed comparably, with perhaps a slight advantage CONCLUSIONS to the quantile normalization. The non-linear method did We have presented three complete data methods of nor- poorer for the spike-in regressions. The scaling method malization and compared these to two different methods had slightly higher slopes. Even so, they were more that make use of a baseline array. Using two different variable. datasets, we established that all three of the complete We saw that the choice of a baseline does have data methods reduced the variation of a probeset measure ramifications on down-stream analysis. Choosing a poor across a set of arrays to a greater degree than the scaling baseline would conceivably give poorer results. We also method and unnormalized data. The non-linear method saw that the complete data methods perform well at both seemed to perform at a level similar to the complete variance reduction and on the matter of bias, and in data methods. Our three complete data methods, while addition more fully reflect the complete set of data. For different, performed comparably at reducing variability this reason we favor a complete data method. across arrays. In terms of speed, for the three complete data methods, When making pairwise comparisons the quantile the quantile method is the fastest. The contrast method method gave the smallest distance between arrays. These is slower and the cyclic loess method is the most time distances also remained fairly constant across intensities. consuming. The contrast and cyclic loess algorithms are In relation to bias, all three complete data methods modifications of an accepted method of normalization. –6 –4 –2 0 2 4 6 Normalization of Oligonucleotide Arrays Hartemink,A., Gifford,D, Jaakkola,T. and Young,R. (2001) Maxi- The quantile method has performed favorably, both in mum likelihood estimation of optimal scaling factors for expres- terms of speed and when using our variance and bias sion array normalization. In SPIE BIOS 2001. criteria, and therefore should be used in preference to the Ihaka,R. and Gentleman,R. (1996) R: a language for data analysis other methods. and graphics. J. Comput. Graph. Stat., 5(3), 299–314. While there might be some advantages to using a Irizarry,R.A., Hobbs,B., Colin, F., Beazer-Barclay,Y.D., Antonellis,K., common, non-data driven, distribution with the quantile Scherf,U. and Speed,T.P. (2003) Exploration, normalization and method, it seems unlikely an agreed standard could be summaries of high density oligonucleotide array probe level data. reached. Different choices of a standard distribution might Biostatistics.in press. be reflected in different estimated fold changes. For this Li,C. and Wong,W.H. (2001) Model-based analysis of oligonu- reason we prefer the minimalist approach of a data based cleotide arrays: model validation, design issues and standard er- normalization. ror applications. Genome Biol., 2(8), 1–11. Lipshutz,R., Fodor,S., Gingeras,T. and Lockart,D. (1999) High den- sity synthetic olignonucleotide arrays. Nature Genet., 21(Suppl), REFERENCES 20–24. Affymetrix (2001) Statistical algorithms reference guide, Technical Schadt,E., Li,C., Su,C. and Wong,W.H. (2001) Analyzing high- report, Affymetrix. density oligonucleotide gene expression array data. J. Cell. Astrand,M. (2001) Normalizing oligonucleotide arrays. Unpub- Biochem., 80, 192–202. lished Manuscript. http://www.math.chalmers.se/∼magnusaa/ Schadt,E., Li,C., Eliss,B. and Wong,W.H. (2002) Feature extraction maffy.pdf and normalization algorithms for high-density oligonucleotide Bolstad,B. (2001) Probe level quantile normalization of high gene expression array data. J. Cell. Biochem., 84(S37), 120–125. density oligonucleotide array data. Unpublished Manuscript. Sidorov,I.A., Hosack,D.A., Gee,D., Yang,J., Cam,M.C., http://www.stat.berkeley.edu/∼bolstad/ Lempicki,R.A. and Dimitrov,D.S. (2002) Oligonucleotide Cleveland,W.S. and Devlin,S.J. (1998) Locally-weighted regres- microarray data distribution and normalization. Information sion: an approach to regression analysis by local fitting. J. Am. Sciences, 146, 65–71. Stat. Assoc., 83, 596–610. Venables,W. and Ripley,B.D. (1997) Modern Applied Statistics with Dudoit,S., Yang,Y.H., Callow,M.J. and Speed,T.P. (2002) Statistical S-PLUS, Second edn, Springer, New York. methods for identifying genes with differential expression in Warrington,J.A., Dee,S. and Trulson,M. (2000) Large-scale replicated cDNA microarray experiments. Stat. Sin., 12(1), 111– genomic analysis using Affymetrix GeneChip .In Schena,M. (ed.), Microarray Biochip Technology. BioTechniques Books, GeneLogic (2002) Datasets http://www.genelogic.com. New York, Chapter 6, pp. 119–148.
Bioinformatics – Oxford University Press
Published: Jan 22, 2003
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.