Bioinformatics Advance Access originally published online on January 17, 2006
Bioinformatics 2006 22(7):802-807; doi:10.1093/bioinformatics/btk047
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Combining signals from spotted cDNA microarrays obtained at different scanning intensities
1Bioinformatics Unit, Institute for Crop Production and Grassland Research, University of Hohenheim Fruwirthstrasse 23, 70599 Stuttgart, Germany
2Center for Plant Molecular Biology, Department of General Genetics, University of Tuebingen Auf der Morgenstelle 28, 72076 Tuebingen, Germany
*To whom correspondence should be addressed.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSALGORITHMIMPLEMENTATIONRESULTSDISCUSSIONREFERENCES |
|---|
Motivation: The analysis of spotted cDNA microarrays involves scanning of color signals from fluorescent dyes. A common problem is that a given scanning intensity is not usually optimal for all spotted cDNAs. Specifically, some spots may be at the saturation limit, resulting in poor separation of signals from different tissues or conditions. The problem may be addressed by multiple scans with varying scanning intensities. Multiple scanning intensities raise the question of how to combine different signals from the same spot, particularly when measurement error is not negligible.
Results: This paper suggests a non-linear latent regression model for this purpose. It corrects for biases caused by the saturation limit and efficiently combines data from multiple scans. Combining multiple scans also allows reduction of technical error particularly for cDNA spots with low signal. The procedure is exemplified using cDNA expression data from maize.
Availability: All methods were implemented using standard procedures available in the SAS/STAT module of the SAS System. Programming statements are available from the first author upon request.
Contact: piepho{at}uni-hohenheim.de
Supplementary information: The supplementary data are available at Bioinformatics online.
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSALGORITHMIMPLEMENTATIONRESULTSDISCUSSIONREFERENCES |
|---|
When scanning fluorescent spot signals (Cy3/Cy5) from cDNA microarrays, the choice of an optimal scanning intensity is critical. A particular choice of scanning intensity may be optimal or near-optimal only for a subset of the spotted cDNAs, while other cDNAs are at the detection limit or saturation limit (Lyng et al., 2004). For this reason, it has been suggested to scan the same microarray repeatedly using different scanning intensities. If there are G genes and S scanning intensities, the procedure results in a G x S array of signals. One possible route for analysis is to select an optimal scanning intensity for each cDNA spot and then base the analysis for a spot only on the signal detected at this one intensity. This approach has several drawbacks. Choice of the optimal scanning intensity is difficult, and some kind of automated analysis is needed. Also, looking at several scanning intensities simultaneously exacerbates the multplicity problem (Korn et al., 2004). Moreover, using only data from one intensity is a waste of information, particularly when error variance is large. A more promising approach is to integrate the results of multiple scans into a single expression value per gene and dye (Romualdi et al., 2003; Garcia de la Nava et al., 2004). An obvious procedure is to compute, for each cDNA spot, some kind of average across scanning intensities. The simple average is not necessarily the best method for doing this, since signals may be more informative for some scanning intensities than for others. Specifically, for spots with signals near the saturation threshold, the information content is relatively low, while signals with intermediate intensity are more informative.
The central idea of this paper is to devise an adjustment for observations near the saturation limit based on a statistical model. The main assumption is that the observed signal at a given scanning intensity and the amount of expression product are connected by some non-linear relationship. The amount of the gene product enters the model as a latent variable, thus giving rise to a non-linear latent regression model. The approach is exemplified using cDNA gene expression data from maize.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSALGORITHMIMPLEMENTATIONRESULTSDISCUSSIONREFERENCES |
|---|
The data
The data are from a large project aimed at unraveling the molecular causes of heterosis in maize. Details of the experimental design are described in Keller et al. (2005). Kernels of the maize inbred lines UH002, UH005, UH250, UH301, and their eight reciprocal hybrids were germinated on filter paper (20 x 70 cm Grade 603 N, Sartorius, Göttingen, Germany) in a phytochamber at 26°C, with a 16 h light, 8 h dark cycle and 60% humidity. RNA was extracted from primary roots 3.5 days after germination using Trizol reagent (Invitrogen, Karlsruhe, Germany) and the Oligotex mRNA Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. mRNA concentrations were quantified with RiboGreen (Molecular Probes, Leiden, Netherlands) using a fluorescence spectrophotometer F-2000 (Hitachi, Japan). mRNA was transcribed into cDNA using random hexamer primer and reverse transcriptase ‘SuperscriptII’ (both Invitrogen). cDNA was indirectly labelled with Cy3 or Cy5 fluorescent dye (Amersham Biosciences, Little Chalfont, UK) and hybridized with a 12.160 element microarray chip (Generation II VersionB) generated at the Iowa State Microarray Facility (Ames, IA; http://www.plantgenomics.iastate.edu/maizechip/), as described in Hedge et al. (2000). We analyzed a total of 49 arrays. The four parents were replicated more often than the hybrids because this optimized accuracy for estimating heterosis (Keller et al., 2005).
Scanning of the cDNA arrays was performed with a GMS 418 Array Scanner (Genetic MicroSystems, Woburn, MA). Each scan was obtained with a laser power of 90%. A series of six scans, in ascending order of PMT gain, was performed at 10 µm resolution. Initial settings for Cy3 and Cy5 were 
10–20% and were increased by 5 or 10% increments. Fluorescent signal intensities were determined and quantified using ImaGene 5.6 (Biodiscovery, Inc., Mariana Del Rey, CA). In total 11 446 spots of this maize microarray chip were subjected to data analysis.
Scanning intensities were log2-transformed. For spots with signals at all six scanning intensities, we computed the mean signal on the log-scale and then plotted each log-signal versus its corresponding spot mean (Fig. 1). The plot suggests that there is a marked saturation level between 15 and 16 units, which is reached for a number of spots at the highest scanning intensities. Also, in the region of linear response, the variance tends to increase towards lower signal intensities. The relationship between any two intensities is linear except for the area around the detection limit. Below the detection limit, the scatter of spots runs parallel for all intensities, and the slope of each scatter is very close to unity.
|
The model
The strong linearity below the saturation limit suggests that up to a region close to the saturation limit it is reasonable to assume a linear relationship between expected signal (µ) and a latent variable (
) related to the amount of gene expression product and the scanning intensity. The parallelism of scatters for the six intensities in Figure 1 further suggests that an additive model of the following form may be used:
 | (1) |
ij is the latent value for the i-th spot at the j-th scanning intensity, 
j is the main effect of the j-th intensity (j = 1, ... , J = 6), subject to the identifiability constraint 6 = 0, and gi is the main effect of the i-th spot (i = 1, ... , I = 11 446). The rationale for this model choice is that a plot of ij versus the spot mean (i1 + i2 + i3 + i4 + i5 + i6)/6 will be linear with slope equal to unity, which corresponds to the pattern in Figure 1.
It appears to be difficult to include the deviant behavior near the saturation point in a single functional model, so one may resort to a segmented regression model (Seber and Wild, 1989; Schabenberger and Pierce, 2002). The simplest such model, also known as the ‘broken stick’ model, has two linear regression lines with different slopes connected at the change point. This has the undesirable property of an abrupt change in slope at the change point. A model with a smooth non-linear transition between linear part and an asymptotic plateau may be obtained by fitting an exponential function near the saturation limit. The model can be written as
 | (2a) |
 | (2b) |
is the saturation limit, 
is the break point, and ß and 
are regression parameters. The model is depicted in Figure 2. The observed signal is modelled as
 | (3) |
|
We require a smooth transition at the break point
ij = in the sense that function value and first derivative of the connected segments must coincide (Seber and Wild, 1989). This requirement gives rise to two constraints on the parameters. Specifically, in order to have a smooth transition in the first derivative, we require
 | (4) |
 | (5) |
j and spot effects gi, only the effects and need to be estimated. |
ALGORITHM |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSALGORITHMIMPLEMENTATIONRESULTSDISCUSSIONREFERENCES |
|---|
One might consider fitting model (2) directly in a single step by non-linear regression. This would be computationally very expensive, however, because the dimension of the problem nearly equals the number of spots (
12 000 in the case of the maize cDNA microarrays used in this paper). The space required to hold and manipulate the design matrix and the data usually exceeds the memory available on standard PCs. Even if memory permits fitting of the model by standard non-linear model procedures, computing time will typically be unacceptably long for routine use. We therefore sought a simple alternative.
The main idea is to fit the model by alternating least squares (Kroonenberg and De Leeuw, 1980), keeping one set of parameters fixed, while estimating the other, and vice versa. The first set of parameters, denoted as intensity set, consists of the intensity effects j and the two non-linear parameters and , while the other set, denoted as spot set, contains the spot effects gi. When the spot set is kept fixed, optimization is with respect to only J + 1 parameters, which is computationally inexpensive. When the intensity set is kept fixed, there are I parameters to be estimated. The key idea is that each of these parameters can be fitted individually, because conditionally on the intensity set, all information on gi is confined to the data of the i-th spot. The algorithm suggested in this paper alternates between the two sets, fitting parameters of the one set, while keeping parameters in the other set fixed. Each of the alternating optimization steps is itself trivial computationally. The total computational cost is governed by the total number of alternating steps. Our experience is that several dozens of iterations are needed until convergence.
The signal versus spot mean plot (Fig. 1) shows some heterogeneity of variance, the variance increasing with decreasing expected signal. Inspection of Figure 1 reveals that heterogeneity is mainly a problem for spots with weak mean signal. One can model the variance as a function of the spot effect and re-estimate all parameters by weighted least squares (WLS), possibly iterating between estimation of the expected signal and the variance (Carroll and Ruppert, 1988). We suggest to use only one or two iterations, following the suggestion by Carroll and Ruppert (1988). The scatterplot of variance (log-scale) versus spot effect may be somewhat irregular (Fig. 3), so we use locally weighted smoothing splines (LOESS; Cleveland, 1979) for fitting the variance function.
|
A further aspect of multiple scan data is that occasionally a signal may be missing. Thus, any proposed method should be able to handle unbalanced data in an efficient way. Standard procedures for linear and non-linear regression (Searle, 1987; Seber and Wild, 1989) are well suited for this task. We initialize spot and intensity effects based on a fit of the additive model in Equation (1). To speed up computing time, intensity effects are estimated by absorbing spot effects into the least squares equations (Engel, 1990; Harville, 1997; SAS Institute, 1999; see Appendix A in the online Supplementary Data). The parameter
, which appears in the exponent of the non-linear response function, is the most sensitive parameter in iterations. Therefore, we determine convergence based on the relative change in this parameter. The following algorithm is proposed to fit model (2):
- Initialize intensity and spot effects by fitting an additive two-way model [Equation (1)] to the data. When fitting the model, absorb spot effects to save computing time. Initialize by max(yij). For , try several starting values between 0 and 5. Initialize weights wi = 1. Set Z = 0. Set Zmax equal to the desired number of WLS steps.
- Keep the spot set (gi) fixed. Estimate the intensity set (, , j) by weighted non-linear least squares using weights wi. Use the whole data set in this optimization.
- Keep the intensity set (, , j) fixed and estimate the spot set (gi) by nonlinear least squares. Perform the optimization separately for each spot.
- If the relative change in the estimate of is smaller than some small value
, proceed to step (5). Otherwise return to step (2).
- If Z = Zmax, terminate the algorithm. Otherwise, do the following. Based on last fit, compute least squares residuals. From the residuals, compute the sample variance per spot. Fit a LOESS curve to predict log-variance from spot effects. Based on the predicted variance, compute weights wi per spot as the inverse of the variance. Set Z
Z + 1 and go back to step (2).
In order to evaluate the algorithm, we compared our method with a standard analysis based on a single scan. For each array, we selected the scan intensity that yielded the best separation among genotypes based on visual inspection, thus mimicking a procedure likely to be applied in practice when multiple scans are available. The data were analyzed using a mixed model approach described in detail by Keller et al. (2005). We tested pairwise contrasts using a Wald-type t-test based on the mixed model fit. We compared P-values of these tests obtained for single-scan data and multiple scans. The P-value, which assesses the significance of a test, is a measure that integrates accuracy of a method: If compared with a standard method P-values of a new method are decreased, this indicates increased accuracy of the new method.
|
IMPLEMENTATION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSALGORITHMIMPLEMENTATIONRESULTSDISCUSSIONREFERENCES |
|---|
We used procedures of the SAS system for all computations (SAS Institute, 1999). Non-linear regression was performed with the NLIN procedure. The LOESS fit for the logarithmically transformed sample variance was computed using the LOESS procedure. Programming statements are available from the first author upon request.
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSALGORITHMIMPLEMENTATIONRESULTSDISCUSSIONREFERENCES |
|---|
The plot of logarithmically transformed variance of residuals versus spot effect shows marked heterogeneity of variance. The plot for example shown in Figure 3 also indicates that the variance stabilizes for high signal values. The weighted estimates of the intensity effects for one array are reported in Table 1. The small standard errors support our claim that intensity effect estimates can approximately be regarded as known. The threshold is estimated as
, while 
. The resulting estimate of the breakpoint is 
. For spots near the saturation limit, we computed estimates 
based on the fitted model and compared these with the simple spot means 
(Fig. 4). It is seen that the model-based adjustment results in an upward shift compared with the naïve spot means. We also plotted estimates of µij versus 
for all six intensities (Fig. 5). The resulting curves closely resemble the plot of yij versus (Fig. 1), suggesting that the model fit is very good.
|
|
|
For 12 contrasts we evaluated the P-values of t-tests based on single and multiple scans. In most cases, multiple scans led to more significances (Table 2). For two contrasts, we plotted P-values of all spots for both analyses. When both tests yield identical result, all spots should fall on a straight line with slope equal to unity. The examples in Figures 6 and 7 show two cases, where P-values of multiple scans were reduced compared with single scans, thus yielding more significant results.
|
|
|
|
DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSALGORITHMIMPLEMENTATIONRESULTSDISCUSSIONREFERENCES |
|---|
The motivation for developing the method proposed in this paper was to address the problem of scanning spots near the saturation limit in the presence of measurement error. Some alternative suggestions have been made on how to use multiple scans to reduce variability of microarray expression data. Our approach is novel in that non-linearity near the saturation limit is explicitly modeled, thus allowing an unbiased adjustment for saturated spots. Also, heterogeneity of variance, which was dramatic in the example used in this paper, is accounted for by a WLS approach. Essentially, our non-linear model linearly extrapolates the observed signal beyond the saturation limit and estimates what would have been observed in the absence of a saturation limit (Garcia de la Nava et al., 2004). The procedure requires multiple scans at different scanning intensities.
A main issue addressed by our method is non-linearity near the saturation limit and the difficulty of determining whether or not a spot signal is at the saturation limit. The difficulty arises from the presence of measurement error. If one could unequivocally identify observations in the saturation region (ij > ), one could delete these and fit the additive model (2a) to the remaining data. Provided for each spot at least one observation obeys the inequality ij > , one could thus obtain unbiased estimates of spot effects gi. In practice, this procedure may be difficult to implement when measurement error is not negligible at the saturation limit, as was the case for the maize data (Fig. 1). Correct classification of an observation as being at the saturation limit is difficult, since both ij and are not directly observable, and one is likely to miss some signals for which ij > , while some signals are kept for which ij < . False classification of observations owing to measurement error may yield a downward bias for spots close to the saturation limit.
Romualdi et al. (2003) propose to use either the maximum or the mean pixel-intensity of several scans for further analysis, arguing that the detection of differentially expressed genes can be considerably enhanced. However, no correction for saturated pixels is provided. Lyng et al. (2004) estimate intensity effects as the mean ratio of untransformed signals at two scanning intensities, assuming a constant ratio across spots not at the saturation limit. This is essentially equivalent to our assumption of additive intensity effects on a logarithmic scale. The approach by Lyng et al. (2004) allows only two scanning intensities, which may not be enough to accommodate optimal scanning intensities for all cDNA spots. Moreover, their mean ratio estimator does not allow to handle heterogeneity of variance (see Appendix B in the online Supplementary Data), which appears to be similarly pronounced in their data as in our maize data. Taking logarithms and working on an additive scale does not fully resolve the problem. We have demonstrated, however, how a simple WLS approach can be implemented to account for heterogeneity of variance.
A correction method for saturation at pixel-level, where a single scanning of an array is sufficient, is provided by Dodd et al. (2004). Similar to Lyng et al. (2004), a linear relationship is estimated, though not between the spots of two scans, but between the pixels of the Cy3 and the Cy5 channel of a certain spot. Thus the pixel of the channel, where no saturation is present, is used for calculating the corrected value for the saturated pixel. The threshold for pixels regarded as saturated is set at a fixed value, which might introduce misclassification, as discussed above. Also, when the number of saturated pixels is high or saturation occurs in both channels, the performance of the procedure might not be optimal.
Garcia de la Nava et al. (2004) provide an approach that adjusts signals not only for saturation but for quantization, i.e. the process of assigning each possible value for the effective amount of hybridization product a discrete value of signal intensity. Each array is scanned at two different intensities. Signal intensities of both scans are related to one another by a linear or a gamma curve using robust regression. With the parameter estimates a maximum likelihood approach is used to gain corrected signal values. Variance heterogeneity does not seem to be an important issue in their data and so is neglected.
This paper has focussed on non-linearities near the saturation limit. In principle, non-linearities could also occur near the detection limit, though this was not a serious issue with our data. Extension of our algorithm to cater for non-linearities near the detection limit is straightforward. One would simply need to add a non-linear segment at low intensities, yielding a model with three segments: the middle segment is linear, while the two outer segments are non-linear. This model would have two additional non-linear parameters compared with model (2).
With multiple laser scans, photobleaching is a potential problem, though we did not observe this effect in our experiments. Photobleaching occurs mainly when scanning intensities are very high, reducing the emitted signal remaining for subsequent scans. As a precaution, we therefore always started with the lowest scanning intensity and finished with the highest. While we did not find evidence of photobleaching with our scans, it should be stressed that our proposed model will automatically account for such effects. Assuming that photobleaching affects all spots in a similar way, the intensity effects j will capture any signal reductions across spots.
Analysis of the maize data revealed that variance increases dramatically as signal intensity decreases (Figs 1 and 3). Thus, repeated scans are also useful for spots with signals well removed from the saturation limit in that they provide measurement replication and thus reduce the standard error of mean signals. Multiple scans are therefore of benefit not only for spots with very high signal, for which saturation is a problem, but also for spots with very low signal, for which measurement error is the primary issue.
Conflict of Interest: none declared.
|
FOOTNOTES |
|---|
Associate Editor: Steen Knudsen
Received on December 9, 2005; revised on December 22, 2005; accepted on January 6, 2006
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSALGORITHMIMPLEMENTATIONRESULTSDISCUSSIONREFERENCES |
|---|
Carroll, R.J. and Ruppert, D. Transformation and Weighting in Regression, (1988) , London Chapman and Hall.
Cleveland, W.S. (1979) Robust locally weighted regression and smoothing scatterplots. J. Am. Stat. Assoc, . 74, 829–836[CrossRef][Web of Science].
Dodd, L.E., et al. (2004) Correcting log ratios for signal saturation in cDNA microarrays. Bioinformatics, 20, 2685–2693
Engel, B. (1990) The analysis of unbalanced linear models with variance components. Stat. Neerl, . 44, 195–219.
Garcia de la Nava, J., et al. (2004) Saturation and quantization reduction in microarray experiments using two scans at different sensitivities. Stat. Appl. Genet. Mol. Biol, . 3, 11.
Harville, D.A. Matrix Algebra from a Statistician's Perspective, (1997) , Berlin Springer.
Hedge, P., et al. (2000) A concise guide to cDNA microarray analysis. Biotechniques, 29, 548–562[Web of Science][Medline].
Keller, B., et al. (2005) Designing a microarray experiment to estimate heterosis. Theoret. Appl. Genet, . 111, 57–64.
Korn, E.L., et al. (2004) Controlling the number of false discoveries: application to high-dimensional genomic data. J. Stat. Plan. Inference, 124, 379–398[CrossRef].
Kroonenberg, P.M. and De Leeuw, J. (1980) Principal component analysis of three-mode data by means of alternating least squares algorithms. Psychometrika, 45, 69–97[CrossRef].
Lyng, H., et al. (2004) Profound influence of microarray scanner characteristics on gene expression ratios: analysis and procedure for correction. BMC Genomics, 5, 10[CrossRef][Medline].
Mood, A., Graybill, F.A., Boes, D.C. Introduction to the Theory of Statistics, (1974) 3rd edn , New York McGraw-Hill.
Romualdi, C., et al. (2003) Improved detection of differentially expressed genes in microarray experiments through multiple scanning and image integration. Nucleic Acids Res, . 31, e149
SAS Institute, Inc. SAS/STAT User's Guide, (1999) , Cary, NC Version 8 SAS Institute.
Schabenberger, O. and Pierce, F.J. Contemporary Statistical Models for the Plant and Soil Sciences, (2002) , Boca Raton CRC Press.
Searle, S.R. Linear Models for Unbalanced Data, (1987) , New York Wiley.
Seber, G.A.F. and Wild, C.J. Nonlinear Regression, (1989) , New York Wiley.
CiteULike 
Connotea 
Del.icio.us What's this?
This article has been cited by other articles:
|
|
 |
|
  N. Hoecker, B. Keller, N. Muthreich, D. Chollet, P. Descombes, H.-P. Piepho, and F. Hochholdinger Comparison of Maize (Zea mays L.) F1-Hybrid and Parental Inbred Line Primary Root Transcriptomes Suggests Organ-Specific Patterns of Nonadditive Gene Expression and Conserved Expression Trends Genetics, July 1, 2008; 179(3): 1275 - 1283. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||