BACKGROUND
microarray technology provides simultaneous measurements of expression levels for thousands of genes.
each step from sample preparation to data analysis, however, contains potential sources of bias and variability.
proper normalization adjusts for differences which interfere with the comparison of intensities of different labels at a given probe and with the comparison of intensities of corresponding probes on different arrays.
proper data normalization should allow for the comparison of expression levels across different arrays.
subsequent data analysis results are heavily dependent on effective normalization.
normalization issues differ for dual-labeled platforms compared to single labeled platforms such as the affymetrix genechip arrays.
in this paper we address normalization for dual-labeled arrays with either cdna or oligonucleotide probes.
the objective of normalization for dual-labeled arrays is to correct for differences in intensities for the two labels on the same array.
these differences arise from factors such as differences in sample concentrations, differences in photomultiplier tube setting, and differences in the affinity of the two labels for dna.
median or mean based global normalization methods use a single normalization factor applied to all genes on the array to adjust for labeling bias  <cit> .
such methods are widely used because of their simplicity.
intensity-based and location-based methods take into account intensity and spatial dependence on dye bias normalization factors  <cit> .
both global and intensity/location based normalization methods assume that most of the genes are not differentially expressed between the two samples hybridized on the array, and that for the differentially expressed genes, the direction of the difference is symmetric between the two samples.
in many important cases, however, these assumptions are not appropriate because: 1) more than half of the genes are differentially expressed on the array; 2) the numbers of over- and under-expressed genes on the array are unequal; 3) only genes of specific biological interest are selected to make a customized array, which are highly variable across the samples.
in the above cases, the global normalization methods and intensity/location based normalization methods become less accurate and a more sophisticated method is needed  <cit> .
there are some methods which attempt to adaptively identify the subset of 'housekeeping' genes  <cit> .
these methods require multiple arrays in order to identify the 'housekeeping' gene set, which does not always exist.
newton et al.
proposed a gamma-gamma-bernoulli model for identifying differentially expressed genes in dual labeled arrays  <cit> .
we have generalized newton's model and here propose an adaptive method based on three-component mixture model for normalization of dual labeled microarray data.
RESULTS
as described in the methods section, we have applied our adaptive method to both the simulated data and real data.
we have also compared our method with the global method and the intensity-based lowess method.
results of the simulation studies are shown as bar plots in figure  <dig>  figure 1a shows the comparison of our adaptive method, the global method and the lowess method when no noise was added.
when the majority of genes in the array were non-differentially expressed , or the numbers of over- and under-expressed genes on the array were equal , the root mean squared error  of the adaptive method was comparable with the other two methods; all were very small.
when the array contained unequal numbers of over- and under-expressed genes and when the majority of genes were differentially expressed , the rmses of the global normalization method and the lowess method were much larger than those of the adaptive method.
the differences ranged from around a two fold difference  when the number of under-, null, and over-expressed genes were  <dig>   <dig>  and  <dig>  to more than a three fold difference  when the number of under-, null, and over-expressed genes were  <dig>   <dig>  and  <dig>  the rmses for the adaptive method ranged from  <dig>  to  <dig>  in log <dig> scale.
we compared the histogram of observed intensities to the fitted marginal density from the adaptive method as a simple check to see whether the proposed model and the estimation procedure are in line with available data.
figure  <dig> shows the histograms of log and log intensities of red and green channels of the simulated data, and the curve in each plot is the estimated density obtained from the fitted model.
it is seen the data fits to the model quite well.
gaussian noise with sd of  <dig>  and  <dig>  were added so that the data was not generated from the same model used for analysis with the adaptive method.
the rmses of the global normalization method and the lowess method remained large, while the rmses of the adaptive method remained small, ranging from  <dig>  to  <dig>  on the log <dig> scale .
in the above simulation, no apparent groups could be seen in the histograms of log .
better results for the adaptive method were also obtained for a simulation case where the three groups  are apparent in the histogram of log.
the results can be seen in figure  <dig> and figure  <dig> .
results comparing rmses for the adaptive method, the global method and the lowess method with real data are shown in table  <dig>  the rmses of the adaptive method on data generated from ten different arrays ranged from  <dig>  to  <dig> , in comparison with rmses of around  <dig>  using the global normalization method.
the average rmse  of the lowess method is almost two times that of our adaptive method , although the lowess method performed better than the global method .
figure  <dig> shows the histograms of log and log intensities of red and green channels of the real data, and the curve in each plot is the estimated density from the adaptive method.
discussion
in this paper, we propose a new method for normalization of dual-labeled arrays in cases where the number of differentially expressed genes is substantial and not necessarily symmetric in direction.
the method performed effectively with both simulated and real data.
we started our model building initially by introducing an unknown constant c into newton's gamma-gamma-bernoulli model  <cit> .
the mixture model consisted of two groups: non-differentially expressed genes  and differentially expressed genes .
log ~ gamma
log ~ gamma
sk ~ gamma
log ~ gamma
log ~ gamma
~ gamma
~ gamma
we found that when the differential expression was symmetric between the two samples, the model worked well.
however, the error increased significantly when the ratio of the numbers of under- to over-expressed genes shifted from  <dig>
in order to make the model more flexible, we modified the model by assigning different scale factors γr and γg for the red channel and green channel intensities.
for this modified two-component mixture, the error still remained large.
we then extended the model into a three-component mixture model listed as equations 8a-8c in the additional material .
the model was then quite flexible but there were too many parameters that needed to be optimized.
after we tested it with simulation data and real data, we found the estimated model was not stable and difficult to optimize.
we finally simplified the model to our final model given by equations 2a-2c .
when applying it to real data or simulated data, the estimates converged well close to globe optima.
when different start points were used, the optimizations remained relatively robust.
evaluation of normalization methods can be difficult since the true normalization factors are unknown with real data for custom arrays.
we avoided this problem by synthesizing customized arrays based on real data for standard arrays containing thousands of genes.
in order to make the distribution of each component group look smoother, we allowed certain range of overlap between the adjacent groups.
additional sampling method was tried to divide the whole distribution range into many non-overlapping intervals.
in each interval the number of genes sampled increased when the absolute value of log <dig> became larger .
the model fitting results using data generated by this sampling method are listed in table  <dig>  and figure  <dig> .
we compared our adaptive method with the global method and the intensity-based lowess method.
the lowess method assumes that in each intensity interval either the majority of genes are non-differentially expressed or the numbers of up- and down-regulated genes are equal.
the global median normalization makes these assumptions only over the array as a whole.
it is not surprising that our method performed much better that the above two methods, because the global median method only works well when the assumptions are valid while the intensity-based lowess method is only effective when there are intensity-dependent biases.
correlation structure is complicated for the thousands of genes on a microarray.
in our model, the intensity of each channel is conditionally independent given the scale parameter, but not marginally independent.therefore, we are not assuming the intensities in two channels are independent.
although we did not generate correlated genes in our simulated data sets, correlations of genes do exist in the real data sets we tested.
spatial correlations are also possible but our method is not designed for that purpose.
yang et al.
proposed using the lowess normalization separately within each grid on the array  <cit> .
our algorithm could be similarly applied within each grid to control for spatial effects.
limited simulations were performed in this study.
we also tried to use real data to test our method.
since an appropriate data set with known normalization factor was not available, we synthesized such data sets by sub-setting large arrays in which the true normalization factor could be accurately estimated.
in the process of synthesizing such small arrays we had to choose an empirical threshold to stratify the differentially expressed genes and non-differentially expressed genes.
although we do not believe that the superiority shown for our algorithm depends critically on the threshold chosen nor on details of the synthesis, it would be preferable to evaluate the algorithm on real data sets with know normalization factors.
although our method is designed for dual-labeled cdna array, it can be extended to single channel affymetrix chip data.
the most popular normalization method for the affymetrix chip compares each array to a single base line array for probe set summaries.
the assumptions behind the normalization method are that the majority of the genes are non-differentially expressed and the numbers of over- and under-expressed genes are roughly equal; the same assumptions as those for dual-labeled cdna arrays.
we could treat the base line array as the 'reference channel' and the other array as the 'test' channel and apply our algorithm to probe set summaries.
for affymetrix chip data, there are multiple base pairs in a probe set and each probe has an intensity measurement.
several alternative normalization methods of affymetrix arrays utilize the probe level information.
for example, method based on an 'invariant set' proposed by li and wong assumes that a probe of a non-differentially expressed genes in two arrays to have similar ranks and uses an iterative procedure to identify the invariant set which presumably consists of points from non-differentially expressed genes  <cit> .
CONCLUSIONS
our new normalization method does not require that the majority of genes be non-differentially expressed, and doesn't require multiple array replicates, dye swaps, spiked controls, or housekeeping genes.
it appears much more effective than standard methods when the numbers of over- and under-expressed genes are unequal, and the majority of the genes are differentially expressed.
it can be very useful for general microarray platforms when samples with very different expression profile are co-hybridized and for custom arrays where the majority of genes are likely to be differentially expressed.
in both of these settings, standard normalization methods are problematic.

 microarray technology provides simultaneous measurements of expression levels for thousands of genes. each step from sample preparation to data analysis, however, contains potential sources of bias and variability. proper normalization adjusts for differences which interfere with the comparison of intensities of different labels at a given probe and with the comparison of intensities of corresponding probes on different arrays. proper data normalization should allow for the comparison of expression levels across different arrays. subsequent data analysis results are heavily dependent on effective normalization. normalization issues differ for dual-labeled platforms compared to single labeled platforms such as the affymetrix genechip arrays. in this paper we address normalization for dual-labeled arrays with either cdna or oligonucleotide probes. the objective of normalization for dual-labeled arrays is to correct for differences in intensities for the two labels on the same array. these differences arise from factors such as differences in sample concentrations, differences in photomultiplier tube setting, and differences in the affinity of the two labels for dna. median or mean based global normalization methods use a single normalization factor applied to all genes on the array to adjust for labeling bias. such methods are widely used because of their simplicity. intensity-based and location-based methods take into account intensity and spatial dependence on dye bias normalization factors. both global and intensity/location based normalization methods assume that most of the genes are not differentially expressed between the two samples hybridized on the array, and that for the differentially expressed genes, the direction of the difference is symmetric between the two samples. in many important cases, however, these assumptions are not appropriate because: 1) more than half of the genes are differentially expressed on the array; 2) the numbers of over- and under-expressed genes on the array are unequal; 3) only genes of specific biological interest are selected to make a customized array, which are highly variable across the samples. in the above cases, the global normalization methods and intensity/location based normalization methods become less accurate and a more sophisticated method is needed. there are some methods which attempt to adaptively identify the subset of'housekeeping' genes. these methods require multiple arrays in order to identify the'housekeeping' gene set, which does not always exist. newton et al. proposed a gamma-gamma-bernoulli model for identifying differentially expressed genes in dual labeled arrays. we have generalized newton's model and here propose an adaptive method based on three-component mixture model for normalization of dual labeled microarray data.  as described in the methods section, we have applied our adaptive method to both the simulated data and real data. we have also compared our method with the global method and the intensity-based lowess method. results of the simulation studies are shown as bar plots in figure figure 1a shows the comparison of our adaptive method, the global method and the lowess method when no noise was added. when the majority of genes in the array were non-differentially expressed, or the numbers of over- and under-expressed genes on the array were equal, the root mean squared error of the adaptive method was comparable with the other two methods; all were very small. when the array contained unequal numbers of over- and under-expressed genes and when the majority of genes were differentially expressed, the rmses of the global normalization method and the lowess method were much larger than those of the adaptive method. the differences ranged from around a two fold difference when the number of under-, null, and over-expressed genes were and to more than a three fold difference when the number of under-, null, and over-expressed genes were and the rmses for the adaptive method ranged from to in log scale. we compared the histogram of observed intensities to the fitted marginal density from the adaptive method as a simple check to see whether the proposed model and the estimation procedure are in line with available data. figure shows the histograms of log and log intensities of red and green channels of the simulated data, and the curve in each plot is the estimated density obtained from the fitted model. it is seen the data fits to the model quite well. gaussian noise with sd of and were added so that the data was not generated from the same model used for analysis with the adaptive method. the rmses of the global normalization method and the lowess method remained large, while the rmses of the adaptive method remained small, ranging from to on the log scale. in the above simulation, no apparent groups could be seen in the histograms of log. better results for the adaptive method were also obtained for a simulation case where the three groups are apparent in the histogram of log. the results can be seen in figure and figure. results comparing rmses for the adaptive method, the global method and the lowess method with real data are shown in table the rmses of the adaptive method on data generated from ten different arrays ranged from to, in comparison with rmses of around using the global normalization method. the average rmse of the lowess method is almost two times that of our adaptive method, although the lowess method performed better than the global method. figure shows the histograms of log and log intensities of red and green channels of the real data, and the curve in each plot is the estimated density from the adaptive method. discussion in this paper, we propose a new method for normalization of dual-labeled arrays in cases where the number of differentially expressed genes is substantial and not necessarily symmetric in direction. the method performed effectively with both simulated and real data. we started our model building initially by introducing an unknown constant c into newton's gamma-gamma-bernoulli model. the mixture model consisted of two groups: non-differentially expressed genes and differentially expressed genes. log ~ gamma log ~ gamma     sk ~ gamma log ~ gamma log ~ gamma     ~ gamma ~ gamma we found that when the differential expression was symmetric between the two samples, the model worked well. however, the error increased significantly when the ratio of the numbers of under- to over-expressed genes shifted from in order to make the model more flexible, we modified the model by assigning different scale factors γr and γg for the red channel and green channel intensities. for this modified two-component mixture, the error still remained large. we then extended the model into a three-component mixture model listed as equations 8a-8c in the additional material. the model was then quite flexible but there were too many parameters that needed to be optimized. after we tested it with simulation data and real data, we found the estimated model was not stable and difficult to optimize. we finally simplified the model to our final model given by equations 2a-2c. when applying it to real data or simulated data, the estimates converged well close to globe optima. when different start points were used, the optimizations remained relatively robust. evaluation of normalization methods can be difficult since the true normalization factors are unknown with real data for custom arrays. we avoided this problem by synthesizing customized arrays based on real data for standard arrays containing thousands of genes. in order to make the distribution of each component group look smoother, we allowed certain range of overlap between the adjacent groups. additional sampling method was tried to divide the whole distribution range into many non-overlapping intervals. in each interval the number of genes sampled increased when the absolute value of log became larger. the model fitting results using data generated by this sampling method are listed in table and figure. we compared our adaptive method with the global method and the intensity-based lowess method. the lowess method assumes that in each intensity interval either the majority of genes are non-differentially expressed or the numbers of up- and down-regulated genes are equal. the global median normalization makes these assumptions only over the array as a whole. it is not surprising that our method performed much better that the above two methods, because the global median method only works well when the assumptions are valid while the intensity-based lowess method is only effective when there are intensity-dependent biases. correlation structure is complicated for the thousands of genes on a microarray. in our model, the intensity of each channel is conditionally independent given the scale parameter, but not marginally independent.therefore, we are not assuming the intensities in two channels are independent. although we did not generate correlated genes in our simulated data sets, correlations of genes do exist in the real data sets we tested. spatial correlations are also possible but our method is not designed for that purpose. yang et al. proposed using the lowess normalization separately within each grid on the array. our algorithm could be similarly applied within each grid to control for spatial effects. limited simulations were performed in this study. we also tried to use real data to test our method. since an appropriate data set with known normalization factor was not available, we synthesized such data sets by sub-setting large arrays in which the true normalization factor could be accurately estimated. in the process of synthesizing such small arrays we had to choose an empirical threshold to stratify the differentially expressed genes and non-differentially expressed genes. although we do not believe that the superiority shown for our algorithm depends critically on the threshold chosen nor on details of the synthesis, it would be preferable to evaluate the algorithm on real data sets with know normalization factors. although our method is designed for dual-labeled cdna array, it can be extended to single channel affymetrix chip data. the most popular normalization method for the affymetrix chip compares each array to a single base line array for probe set summaries. the assumptions behind the normalization method are that the majority of the genes are non-differentially expressed and the numbers of over- and under-expressed genes are roughly equal; the same assumptions as those for dual-labeled cdna arrays. we could treat the base line array as the'reference channel' and the other array as the'test' channel and apply our algorithm to probe set summaries. for affymetrix chip data, there are multiple base pairs in a probe set and each probe has an intensity measurement. several alternative normalization methods of affymetrix arrays utilize the probe level information. for example, method based on an'invariant set' proposed by li and wong assumes that a probe of a non-differentially expressed genes in two arrays to have similar ranks and uses an iterative procedure to identify the invariant set which presumably consists of points from non-differentially expressed genes.  our new normalization method does not require that the majority of genes be non-differentially expressed, and doesn't require multiple array replicates, dye swaps, spiked controls, or housekeeping genes. it appears much more effective than standard methods when the numbers of over- and under-expressed genes are unequal, and the majority of the genes are differentially expressed. it can be very useful for general microarray platforms when samples with very different expression profile are co-hybridized and for custom arrays where the majority of genes are likely to be differentially expressed. in both of these settings, standard normalization methods are problematic. 
we propose a mixture model based normalization method that adaptively identifies non-differentially expressed genes and thereby substantially improves normalization for dual-labeled arrays in settings where the assumptions of global normalization are problematic. the new normalization method is effective for general microarray platforms when samples with very different expression profile are co-hybridized and for custom arrays where the majority of genes are likely to be differentially expressed. the new method is evaluated using both simulated and real data. these assumptions can be inappropriate for custom arrays or arrays in which the reference rna is very different from the experimental samples.
 normalization is a critical step in analysis of gene expression profiles. for dual-labeled arrays, global normalization assumes that the majority of the genes on the array are non-differentially expressed between the two channels and that the number of over-expressed genes approximately equals the number of under-expressed genes. these assumptions can be inappropriate for custom arrays or arrays in which the reference rna is very different from the experimental samples.  we propose a mixture model based normalization method that adaptively identifies non-differentially expressed genes and thereby substantially improves normalization for dual-labeled arrays in settings where the assumptions of global normalization are problematic. the new method is evaluated using both simulated and real data.  the new normalization method is effective for general microarray platforms when samples with very different expression profile are co-hybridized and for custom arrays where the majority of genes are likely to be differentially expressed. 