the statistical model proposed in   for one probe set in multiple oligonucleotide arrays has the form it states that the perfect match ( pm)/mismatch ( mm ) difference in array i , probe j of this probe set is the product of model - based expression index ( mbei ) in array i ( i ) and probe - sensitivity index of probe j ( j ) plus random error . fitting the model , we can identify cross - hybridizing probes ( j with large standard error ( se ) , which are excluded during iterative fitting ) and arrays with image contamination at this probe set ( i with large se ) , as well as single outliers ( image spikes ) which are replaced by the fitted values . in effect the estimated expression index i is a weighted average of pm / mm differences : with larger weights given to probes with larger . the image of outliers ( array and single outliers ) identified through model - fitting can be used to assess the quality of an experiment and to identify unexpected problems such as a misaligned corner of a dat file . we have investigated several important properties of the model , including the reliability and stability of the fitted parameters mbei (  ) and probe sensitivity indexes ( ) , the performance of mbei compared to the commonly used average difference ( ad ) , and how the availability of se facilitates downstream comparative and clustering analysis . in practice , in an array experiment , a researcher hybridizes tissue or cell line samples , corresponding to different treatments or conditions , to a batch of arrays . ideally , the probe - sensitivity index ( ) should be independent of the tissue type . this condition , however , may not hold for those probes that have cross - hybridization affinity to non - target genes . nevertheless , assuming that a non - target gene cross - hybridizes only to a few probes of a probe set , and its expression levels across arrays do not correlate with the target gene , the iterative probe - excluding procedure in   may be able to exclude cross - hybridizing probes , regardless of the tissue type hybridized . in addition , the relative probe - sensitivity indexes of the good probes called by the model are likely to be similar across sets of arrays hybridizing to different tissue samples . we apply the model ( equation 1 ) independently to six sets of hu6800 arrays ( 21 leukemia , lymphoma and mantle cell samples , 20 prostate cancer cell lines , 17 brain tumor samples , 55 cancer cell lines , 58 brain samples , and 55 lung tumor samples ) . figure 1a shows the   values fitted for probe set 6457 ( used in figure 1 and 2 of ) in the six array sets . the   patterns resemble each other greatly , showing that the probe - sensitivity index is an inherent property of these non - cross - hybridizing probes and can be consistently identified from different sets of arrays . it is noteworthy that the probe 11 in array set 5 is likely to be cross - hybridizing , making its relative strength ( here mm is consistently larger than pm and this leads to a negative ) dissimilar to the probe 11 in other array sets . the model identifies this probe as a ' probe - outlier ' only for array set 5 and excludes it when calculating mbei ( 0 ) for array set 5 . values estimated for probe sets ( a ) 6457 , ( b ) 1248 , and ( c ) 6571 in six array sets ( shown in panels 16 from left to right for each probe set ) .   values ( constrained to have sum square equal to number of probes used in each array set ) are on the y - axis , and probe pairs are labeled 1 to 20 on the x - axis . the title of each panel ( for example , p = 0 ) indicates the proportion of arrays ' present ' for the target gene in the array set . they are stratified by average lower presence proportion in two array sets ( the presence proportion of a probe set is the proportion of arrays in an array set where the target gene is called ' present ' by genechip 's algorithm ) . the average is taken over c(6 , 2 ) = 15 pairwise comparison of two array sets for each probe set , and the correlation is calculated using probes that are not identified as an outlier in both array sets . the range of the average lower presence proportion for the six boxplots are : ( 0 , 0.17 ) , ( 0.17 , 0.34 ) , ( 0.34 , 0.51 ) , ( 0.51 , 0.68 ) , ( 0.68 , 0.85 ) , ( 0.85 , 1 ) . the title of each boxplot is the number of probe sets classified into this boxplot . eleven probe sets with too few non - outlier probes to calculate   correlations for all 15 comparisons are not included in the boxplots . the average lower presence proportion and average pairwise correlation for probe sets in figure 1 are ( a ) 1 , 0.95 ; ( b ) , 0.93 , 0.94 ; and ( c ) 0 , 0.86 . in figure 1a , b the target gene is present in most samples of all array sets . for a probe set whose target gene is mostly absent throughout samples ( figure 1c ) , many probes are identified as probe - outliers because of their negative indexes . here , we can not obtain correct probe - sensitivity indexes because of the absence of the target gene . nevertheless , the pm - mm values for these probes are random fluctuations around zero , leading to a correct expression index close to zero . if the target gene becomes available for a future array set , the correct probe - sensitivity indexes will be recovered and these probes will be used for expression calculation . occasionally , a responsive probe set may give rise to very different   estimates in two array sets . in figure 1b , probes 8 and 13 have different relative responses in array set 1 and 4 , leading to different probe - response patterns . this might be due to the possibility that the probes in this probe set are differentially cross - hybridized in different array sets , or that the same probe in different batches of arrays may systematically behave differently . identification and flagging such probe sets is desirable and essential if we want to compare arrays hybridized to different tissue samples . figure 2 shows the boxplots of average pairwise correlations of   values between two array sets , stratified by average lower presence proportion in the two sets . in general , when a gene is present in many samples of two array sets , the   patterns estimated from the two sets are very similar . the target gene 's presence in many arrays of an array set allows the probe - sensitivity index to be estimated accurately . from figure 1 of , one can see that some mm probes may respond poorly to the changes in the expression level of the target gene . this phenomenon raised questions on the efficiency of using mm probes , and led some investigators to design custom arrays that use pm probes exclusively ( r. abagyan and yingyao zhou , personal communication ; b.r . conklin , personal communication ) , and others to calculate fold changes using only pm probes ( f. naef , personal communication ) . this design greatly increases the number of genes that can be studied on one array . to investigate the relative performance of pm / mm versus pm - only designs , we exploited the model to estimate gene expression levels using only pm probes , and compared it to the mbei using both pm and mm probes . the full intensity model ( equation 1 of ) specifies the relationship of pm probe responses and expression level  : where vj is the baseline response of probe pair j due to nonspecific hybridization , and ' j is the sensitivity of pm probe of the probe pair j. the parameter estimates can be obtained by iteratively fitting i and vj,'j , regarding the other set as known . the mm probe responses have a similar form as equation 2 except for different probe - sensitivity indexes . we fit a pm - only and an mm - only model to obtain expression values of all 20-probe probe sets using array set 1 . for comparison , we also used half of the probe pairs ( by alternatively picking one out of every two probes ) in a 20-probe probe set to fit to the difference model ( equation 1 ) . for each probe set , these three sets of expression values were compared with the expression values of the original difference model using 20 probes , in terms of correlation of s obtained by two methods across the 21 arrays . we assumed the 20-probe difference model provides the most accurate expression estimates . if , for a probe set , a simplified model ( pm - only , mm - only or 10-probe difference model ) performs reasonably well , we expect its  estimates to correlate with that from the 20-probe difference model . figure 3 shows the histogram and figure 4 the boxplot of correlations of s estimated from the 20-probe difference model and s estimated from the 10-probe difference model ( a ) , the 20-probe pm - only model ( b ) and the 20-probe mm - only model ( c ) . for probe sets with high presence proportion , both the 10-probe difference model and the pm - only model correlate well with the 20-probe difference model . we note that this comparison is intrinsically biased in favor of the 10-probe difference model because the ' truth ' is constructed from pm - mm differences . histogram of correlations between model - based expression values estimated using the 20-probe difference model and those estimated using different models . ( a ) 10-probe difference model ; ( b ) 20-probe pm - only model ; ( c ) 20-probe mm - only model . boxplot of correlations between  values estimated using the 20-probe difference model and s estimated using different models , stratified by presence proportion . ( a ) 10-probe difference model ; ( b ) 20-probe pm - only model ; and ( c ) 20-probe mm - only model . the number of presence calls for a probe set in the 21 arrays and the subpopulation size for the six boxplots are : 03 , 4,385 ; 47 , 693 ; 811 , 413 ; 1215 , 488 ; 1619 , 497 ; and 2021 , 323 . this comparison corroborates the basic notion of the technology : the pm probes hybridize more strongly to the target signals than mm probes and contain most of the information . we stress that , whereas the above analysis illustrates the applicability of model - based analysis to pm - only arrays , the assessment presented here is only tentative because of the limited information provided by the hu6800 arrays on the comparisons . definitive comparisons of the efficiency of the designs must await the availability of data from pm - only arrays . each pair consists of two arrays hybridizing to samples replicated at total mrna level ( the total mrna sample is split and then amplified and labeled separately , and hybridized to two different arrays ) . the differences between the expression values of the two replicate arrays in a pair are due to the variation introduced in experimental steps after the split , the array manufacturing difference and analytical methods such as normalization and expression calculation . this difference provides a lower bound of biological variation that can be detected between two independently amplified samples , and serves as a good statistic for comparing different analytical methods . for comparison , we also used the method in   to calculate ads for all probe sets and plot them in figure 5b ( ad is based on normalized probe values , see methods and materials section for the normalization method . also note that genechip software excludes probes whose pm / mm difference is outside three standard deviations ( sds ) of all probe differences in either of the two arrays in the comparison ; here , as we are comparing multiple arrays at the same time , when calculating ads a probe is excluded if its difference is an outlier in the above sense in any of the arrays , until a minimum of five probes is reached , where all five probes will be used ) . both the mbei and the ad method yielded some expression values differing by more than a factor of two , especially for genes at low expression level . this might be explained by the relatively larger amplification variation for weakly expressed genes , given a constant success rate of amplifying a sequence by a certain fold . log ( base 10 ) expression indexes of a pair of replicate arrays ( array 1 and 2 of array set 5 ) for different statistical methods . only 6,695 ( a ) and 4,696 ( b ) probe sets with positive values in both arrays are used . the center line is y = x , and the flanking lines indicate the difference of a factor of two . researchers often use ' log ratio ' between expression values of a gene in two arrays as the criterion for identifying differentially expressed genes . between duplicate arrays , we expect these log ratios of expression values based on a good expression index ( ad or mbei ) to be close to zero . thus for every probe set we calculated its average absolute log ( base 10 ) ratio of 29 pairs of duplicates as a statistic to compare the variation in expression levels between duplicates using the ad or the mbei method . the average absolute log ratio distribution of the mbei method is significantly lower than that of the ad method when expression level is low ( and thus probe sets have a low proportion of detections of the target gene across arrays ) . as expression level becomes higher ( when the target gene of a probe set is detected in more arrays ) , the ad method shows a rapid improvement in performance , approaching the level of the mbei method . the same boxplots ( figure 7 ) for another set of 60 human u95a arrays consisting of 30 replicate pairs conveys similar information . these results suggest that the mbei method is able to extend the reliable detection limit of expression to a lower mrna concentration . boxplots of average absolute log ( base 10 ) ratios between replicate arrays stratified by presence proportion for different statistical methods . ( a ) mbei method ; ( b ) ad method . the number of presence calls for a probe set in the 58 arrays for the six boxplots are : 09 , 1019 , 2029 , 3039 , 4049 , 5058 . the title of each boxplot is the number of probe sets used for the boxplot . log ratios are not calculated for negative expression values or expression values identified as ' array - outliers ' by the mbei method in either array of a replicate pair , and are not used to calculate the average . 744 probe sets are not included as their average absolute log ratios can not be calculated for all the 29 pairs using either method . similar plots as in figure 6 for another set of 30 pairs of duplicated human u95a arrays . ( a ) mbei method ; ( b ) ad method.the number of presence calls for a probe set in the 60 arrays for the six boxplots are : 09 , 1019 , 2029 , 3039 , 4049 , 5060 . the title of each boxplot is the number of probe sets used for the boxplot . after obtaining expression indexes using ad or mbei , fold changes can be calculated between two arrays for every gene and used to identify differentially expressed genes . usually , low or negative expressions are truncated to a small number before calculating fold changes , and genechip also cautions against using fold changes when the baseline expression is absent . the availability of ses for the model - based expression indexes allows us to obtain confidence intervals for fold changes . suppose where 1 and 2 are the real expression levels in the sample , and 1 and 2 are the model - based estimates of expression levels . we substitute the model - based ses for 1 and 2 letting r = 1/2 be the real fold change , then inference on r can be based on the quantity it can be shown that q has a  distribution with 1 degree of freedom irrespective of the values of 1 and 2 . thus q is a pivotal quantity involving r. we can use q to construct fixed - level tests and to invert them to obtain confidence intervals ( ci ) for fold changes . table 1 presents the estimated expression indexes ( with ses ) in two arrays and the 90% confidence intervals of the fold changes for 14 genes . although all genes have similar estimated fold changes , the confidence intervals are very different . for example , gene 1 has fold change 2.47 and a tight confidence interval ( 2.06 , 3.02 ) . in contrast , gene 11 has a similar fold change of 2.48 but a much wider confidence interval ( 0.96 , 18.18 ) . thus the fold change around 2.5 for gene 11 is not as trustworthy as that for gene 1 . further examination reveals that this is due to the large ses relative to the expression indexes for gene 11 . this agrees with the intuition that when one or both expression levels are close to zero for one gene , the fold change can not be estimated with much accuracy . in addition , when image contamination results in unreliable expression values with large ses , the fold changes calculated using these expression value are attached with wide cis . in this manner , the measurement accuracy of expression values propagates to the estimation of fold changes . using expression levels and associated ses to determine confidence intervals of fold changes in practice , we find it useful to sort genes by the lower confidence bound ( ' lower cb ' in table 1 ) , which is a conservative estimate of the fold change . when an expression index is negative ( as a result of taking pm / mm differences ) , we do not calculate the confidence intervals . in such a case , it is more helpful to filter genes by presence calls . cluster analysis is a popular method for analyzing the data of a series of microarrays [ 6 , 7 ] . if two genes are co - regulated at the transcription level , their expression values across samples are likely to be correlated . clustering algorithms use these correlations ( or monotone transformation of correlations ) to cluster co - regulated genes together the correlation based on the estimated expression levels may , however , be different from that based on the real but unobserved expression levels . also , the commonly used hierarchical clustering algorithm is an irreversible process : once two genes or nodes are merged , they will stay together , even if later on there is good reason to adjust previous clustering . a global way of using se in hierarchical clustering is to resample or bootstrap   the whole ' gene by sample ' data matrix and redo the clustering , then investigate the overall properties emerging from this repertoire of clustering trees . in bittner et al . , the data matrix coming from cdna microarray experiments is resampled using the estimated variation derived from the median sd of log ratios for a gene across samples . as we now have ses for all data points , we can resample each expression value from a normal distribution with mean equal to the estimated expression value and sd equal to the attached se . figure 8a shows a hierarchical clustering tree of 225 selected genes with presence proportion > 0.5 and coefficient of variation ( sd / mean ) > 0.7 across the 20 samples in array set 2 . in trying to interpret this tree , we may be interested in the gene cluster colored in blue and the reliability of the gene members belonging to this cluster . the whole data matrix is resampled , and the clustering is performed again ( figure 8b ) . we notice that some blue genes ( genes in the original cluster are colored blue ) are clustered with other non - blue genes , and some non - blue genes are mixed into the main body of the blue genes . after each resampling , we identify a cluster that contains more than 80% of all the blue genes , but as few non - blue genes as possible ( measured as a percentage of all genes in this cluster ) . this cluster is considered to be the cluster that corresponds to the original one in figure 8a . in figure 8b the root node of the ' corresponding cluster ' is marked with small horizontal line intersecting the vertical line ( representing the range of the cluster ) on the right of the clustering picture . then , for each of all the 225 genes , if it belongs to this ' corresponding cluster ' , we increase its ' in - cluster ' count by 1 . after resampling 30 times , the in - cluster counts are indicated in gray - scale on the left side of the original clustering ( figure 8c ) , with black representing 30 and white representing zero . a high ' in - cluster ' count indicates a gene ' remains ' in the original cluster in most of the resampled clustering trees . gene clustering ( a ) 225 filtered genes are clustered based on their expression profiles across 20 samples . each gene 's expression values are standardized to have mean 0 and sd 1 across 20 samples . although the original ' blue ' genes are scattered to various places , we can still determine where the original cluster is , using the criteria described in the text . ( c ) after resampling 30 times , the reliability of the genes belonging to the original cluster is indicated by the vertical gray - scale bar on the left of the blue - red picture . we can see from figure 8c that most genes in the original cluster are reliable members , whereas a few genes at the bottom of the cluster are not ( in fact they are merged into the original cluster last ) . interestingly , some genes originally not in the original cluster group with the ' corresponding clusters ' during resampling many times and have gray ' in - cluster ' marks . these genes may be related to the original cluster in some way . in summary , this method can help us to distinguish reliable and unreliable gene members of a cluster , as well as draw our attention to related genes originally clustered somewhere else because of the accidental nature of hierarchical clustering . in practice , in an array experiment , a researcher hybridizes tissue or cell line samples , corresponding to different treatments or conditions , to a batch of arrays . ideally , the probe - sensitivity index ( ) should be independent of the tissue type . this condition , however , may not hold for those probes that have cross - hybridization affinity to non - target genes . nevertheless , assuming that a non - target gene cross - hybridizes only to a few probes of a probe set , and its expression levels across arrays do not correlate with the target gene , the iterative probe - excluding procedure in   may be able to exclude cross - hybridizing probes , regardless of the tissue type hybridized . in addition , the relative probe - sensitivity indexes of the good probes called by the model are likely to be similar across sets of arrays hybridizing to different tissue samples . we apply the model ( equation 1 ) independently to six sets of hu6800 arrays ( 21 leukemia , lymphoma and mantle cell samples , 20 prostate cancer cell lines , 17 brain tumor samples , 55 cancer cell lines , 58 brain samples , and 55 lung tumor samples ) . figure 1a shows the   values fitted for probe set 6457 ( used in figure 1 and 2 of ) in the six array sets . the   patterns resemble each other greatly , showing that the probe - sensitivity index is an inherent property of these non - cross - hybridizing probes and can be consistently identified from different sets of arrays . it is noteworthy that the probe 11 in array set 5 is likely to be cross - hybridizing , making its relative strength ( here mm is consistently larger than pm and this leads to a negative ) dissimilar to the probe 11 in other array sets . the model identifies this probe as a ' probe - outlier ' only for array set 5 and excludes it when calculating mbei ( 0 ) for array set 5 . values estimated for probe sets ( a ) 6457 , ( b ) 1248 , and ( c ) 6571 in six array sets ( shown in panels 16 from left to right for each probe set ) .   values ( constrained to have sum square equal to number of probes used in each array set ) are on the y - axis , and probe pairs are labeled 1 to 20 on the x - axis . the title of each panel ( for example , p = 0 ) indicates the proportion of arrays ' present ' for the target gene in the array set . they are stratified by average lower presence proportion in two array sets ( the presence proportion of a probe set is the proportion of arrays in an array set where the target gene is called ' present ' by genechip 's algorithm ) . the average is taken over c(6 , 2 ) = 15 pairwise comparison of two array sets for each probe set , and the correlation is calculated using probes that are not identified as an outlier in both array sets . the range of the average lower presence proportion for the six boxplots are : ( 0 , 0.17 ) , ( 0.17 , 0.34 ) , ( 0.34 , 0.51 ) , ( 0.51 , 0.68 ) , ( 0.68 , 0.85 ) , ( 0.85 , 1 ) . the title of each boxplot is the number of probe sets classified into this boxplot . eleven probe sets with too few non - outlier probes to calculate   correlations for all 15 comparisons are not included in the boxplots . the average lower presence proportion and average pairwise correlation for probe sets in figure 1 are ( a ) 1 , 0.95 ; ( b ) , 0.93 , 0.94 ; and ( c ) 0 , 0.86 . in figure 1a , b the target gene is present in most samples of all array sets . for a probe set whose target gene is mostly absent throughout samples ( figure 1c ) , many probes are identified as probe - outliers because of their negative indexes . here , we can not obtain correct probe - sensitivity indexes because of the absence of the target gene . nevertheless , the pm - mm values for these probes are random fluctuations around zero , leading to a correct expression index close to zero . if the target gene becomes available for a future array set , the correct probe - sensitivity indexes will be recovered and these probes will be used for expression calculation . occasionally , a responsive probe set may give rise to very different   estimates in two array sets . in figure 1b , probes 8 and 13 have different relative responses in array set 1 and 4 , leading to different probe - response patterns . this might be due to the possibility that the probes in this probe set are differentially cross - hybridized in different array sets , or that the same probe in different batches of arrays may systematically behave differently . identification and flagging such probe sets is desirable and essential if we want to compare arrays hybridized to different tissue samples . figure 2 shows the boxplots of average pairwise correlations of   values between two array sets , stratified by average lower presence proportion in the two sets . in general , when a gene is present in many samples of two array sets , the   patterns estimated from the two sets are very similar . the target gene 's presence in many arrays of an array set allows the probe - sensitivity index to be estimated accurately . from figure 1 of , one can see that some mm probes may respond poorly to the changes in the expression level of the target gene . this phenomenon raised questions on the efficiency of using mm probes , and led some investigators to design custom arrays that use pm probes exclusively ( r. abagyan and yingyao zhou , personal communication ; b.r . conklin , personal communication ) , and others to calculate fold changes using only pm probes ( f. naef , personal communication ) . this design greatly increases the number of genes that can be studied on one array . to investigate the relative performance of pm / mm versus pm - only designs , we exploited the model to estimate gene expression levels using only pm probes , and compared it to the mbei using both pm and mm probes . the full intensity model ( equation 1 of ) specifies the relationship of pm probe responses and expression level  : where vj is the baseline response of probe pair j due to nonspecific hybridization , and ' j is the sensitivity of pm probe of the probe pair j. the parameter estimates can be obtained by iteratively fitting i and vj,'j , regarding the other set as known . the same outlier exclusion procedure in   is applied . the mm probe responses have a similar form as equation 2 except for different probe - sensitivity indexes . we fit a pm - only and an mm - only model to obtain expression values of all 20-probe probe sets using array set 1 . for comparison , we also used half of the probe pairs ( by alternatively picking one out of every two probes ) in a 20-probe probe set to fit to the difference model ( equation 1 ) . for each probe set , these three sets of expression values were compared with the expression values of the original difference model using 20 probes , in terms of correlation of s obtained by two methods across the 21 arrays . if , for a probe set , a simplified model ( pm - only , mm - only or 10-probe difference model ) performs reasonably well , we expect its  estimates to correlate with that from the 20-probe difference model . figure 3 shows the histogram and figure 4 the boxplot of correlations of s estimated from the 20-probe difference model and s estimated from the 10-probe difference model ( a ) , the 20-probe pm - only model ( b ) and the 20-probe mm - only model ( c ) . for probe sets with high presence proportion , both the 10-probe difference model and the pm - only model correlate well with the 20-probe difference model . the mm - only model yields noticeably lower correlations , however . we note that this comparison is intrinsically biased in favor of the 10-probe difference model because the ' truth ' is constructed from pm - mm differences . histogram of correlations between model - based expression values estimated using the 20-probe difference model and those estimated using different models . ( a ) 10-probe difference model ; ( b ) 20-probe pm - only model ; ( c ) 20-probe mm - only model . boxplot of correlations between  values estimated using the 20-probe difference model and s estimated using different models , stratified by presence proportion . ( a ) 10-probe difference model ; ( b ) 20-probe pm - only model ; and ( c ) 20-probe mm - only model . the number of presence calls for a probe set in the 21 arrays and the subpopulation size for the six boxplots are : 03 , 4,385 ; 47 , 693 ; 811 , 413 ; 1215 , 488 ; 1619 , 497 ; and 2021 , 323 . this comparison corroborates the basic notion of the technology : the pm probes hybridize more strongly to the target signals than mm probes and contain most of the information . we stress that , whereas the above analysis illustrates the applicability of model - based analysis to pm - only arrays , the assessment presented here is only tentative because of the limited information provided by the hu6800 arrays on the comparisons . definitive comparisons of the efficiency of the designs must await the availability of data from pm - only arrays . each pair consists of two arrays hybridizing to samples replicated at total mrna level ( the total mrna sample is split and then amplified and labeled separately , and hybridized to two different arrays ) . the differences between the expression values of the two replicate arrays in a pair are due to the variation introduced in experimental steps after the split , the array manufacturing difference and analytical methods such as normalization and expression calculation . this difference provides a lower bound of biological variation that can be detected between two independently amplified samples , and serves as a good statistic for comparing different analytical methods . for comparison , we also used the method in   to calculate ads for all probe sets and plot them in figure 5b ( ad is based on normalized probe values , see methods and materials section for the normalization method . also note that genechip software excludes probes whose pm / mm difference is outside three standard deviations ( sds ) of all probe differences in either of the two arrays in the comparison ; here , as we are comparing multiple arrays at the same time , when calculating ads a probe is excluded if its difference is an outlier in the above sense in any of the arrays , until a minimum of five probes is reached , where all five probes will be used ) . both the mbei and the ad method yielded some expression values differing by more than a factor of two , especially for genes at low expression level . this might be explained by the relatively larger amplification variation for weakly expressed genes , given a constant success rate of amplifying a sequence by a certain fold . log ( base 10 ) expression indexes of a pair of replicate arrays ( array 1 and 2 of array set 5 ) for different statistical methods . ( a ) mbei method ; ( b ) ad method . only 6,695 ( a ) and 4,696 ( b ) probe sets with positive values in both arrays are used . the center line is y = x , and the flanking lines indicate the difference of a factor of two . researchers often use ' log ratio ' between expression values of a gene in two arrays as the criterion for identifying differentially expressed genes . between duplicate arrays , we expect these log ratios of expression values based on a good expression index ( ad or mbei ) to be close to zero . thus for every probe set we calculated its average absolute log ( base 10 ) ratio of 29 pairs of duplicates as a statistic to compare the variation in expression levels between duplicates using the ad or the mbei method . the average absolute log ratio distribution of the mbei method is significantly lower than that of the ad method when expression level is low ( and thus probe sets have a low proportion of detections of the target gene across arrays ) . as expression level becomes higher ( when the target gene of a probe set is detected in more arrays ) , the ad method shows a rapid improvement in performance , approaching the level of the mbei method . the same boxplots ( figure 7 ) for another set of 60 human u95a arrays consisting of 30 replicate pairs conveys similar information . these results suggest that the mbei method is able to extend the reliable detection limit of expression to a lower mrna concentration . boxplots of average absolute log ( base 10 ) ratios between replicate arrays stratified by presence proportion for different statistical methods . ( a ) mbei method ; ( b ) ad method . the number of presence calls for a probe set in the 58 arrays for the six boxplots are : 09 , 1019 , 2029 , 3039 , 4049 , 5058 . the title of each boxplot is the number of probe sets used for the boxplot . log ratios are not calculated for negative expression values or expression values identified as ' array - outliers ' by the mbei method in either array of a replicate pair , and are not used to calculate the average . 744 probe sets are not included as their average absolute log ratios can not be calculated for all the 29 pairs using either method . similar plots as in figure 6 for another set of 30 pairs of duplicated human u95a arrays . ( a ) mbei method ; ( b ) ad method.the number of presence calls for a probe set in the 60 arrays for the six boxplots are : 09 , 1019 , 2029 , 3039 , 4049 , 5060 . the title of each boxplot is the number of probe sets used for the boxplot . after obtaining expression indexes using ad or mbei , fold changes can be calculated between two arrays for every gene and used to identify differentially expressed genes . usually , low or negative expressions are truncated to a small number before calculating fold changes , and genechip also cautions against using fold changes when the baseline expression is absent . the availability of ses for the model - based expression indexes allows us to obtain confidence intervals for fold changes . suppose where 1 and 2 are the real expression levels in the sample , and 1 and 2 are the model - based estimates of expression levels . we substitute the model - based ses for 1 and 2 letting r = 1/2 be the real fold change , then inference on r can be based on the quantity it can be shown that q has a  distribution with 1 degree of freedom irrespective of the values of 1 and 2 . thus q is a pivotal quantity involving r. we can use q to construct fixed - level tests and to invert them to obtain confidence intervals ( ci ) for fold changes . table 1 presents the estimated expression indexes ( with ses ) in two arrays and the 90% confidence intervals of the fold changes for 14 genes . although all genes have similar estimated fold changes , the confidence intervals are very different . for example , gene 1 has fold change 2.47 and a tight confidence interval ( 2.06 , 3.02 ) . in contrast , gene 11 has a similar fold change of 2.48 but a much wider confidence interval ( 0.96 , 18.18 ) . thus the fold change around 2.5 for gene 11 is not as trustworthy as that for gene 1 . further examination reveals that this is due to the large ses relative to the expression indexes for gene 11 . this agrees with the intuition that when one or both expression levels are close to zero for one gene , the fold change can not be estimated with much accuracy . in addition , when image contamination results in unreliable expression values with large ses , the fold changes calculated using these expression value are attached with wide cis . in this manner , the measurement accuracy of expression values propagates to the estimation of fold changes . using expression levels and associated ses to determine confidence intervals of fold changes in practice , we find it useful to sort genes by the lower confidence bound ( ' lower cb ' in table 1 ) , which is a conservative estimate of the fold change . when an expression index is negative ( as a result of taking pm / mm differences ) , we do not calculate the confidence intervals . in such a case , it is more helpful to filter genes by presence calls . cluster analysis is a popular method for analyzing the data of a series of microarrays [ 6 , 7 ] . if two genes are co - regulated at the transcription level , their expression values across samples are likely to be correlated . clustering algorithms use these correlations ( or monotone transformation of correlations ) to cluster co - regulated genes together . the correlation based on the estimated expression levels may , however , be different from that based on the real but unobserved expression levels . also , the commonly used hierarchical clustering algorithm is an irreversible process : once two genes or nodes are merged , they will stay together , even if later on there is good reason to adjust previous clustering . thus there is a need to assess the reliability of clusters . a global way of using se in hierarchical clustering is to resample or bootstrap   the whole ' gene by sample ' data matrix and redo the clustering in bittner et al . , the data matrix coming from cdna microarray experiments is resampled using the estimated variation derived from the median sd of log ratios for a gene across samples . as we now have ses for all data points , we can resample each expression value from a normal distribution with mean equal to the estimated expression value and sd equal to the attached se . figure 8a shows a hierarchical clustering tree of 225 selected genes with presence proportion > 0.5 and coefficient of variation ( sd / mean ) > 0.7 across the 20 samples in array set 2 . in trying to interpret this tree , we may be interested in the gene cluster colored in blue and the reliability of the gene members belonging to this cluster . the whole data matrix is resampled , and the clustering is performed again ( figure 8b ) . we notice that some blue genes ( genes in the original cluster are colored blue ) are clustered with other non - blue genes , and some non - blue genes are mixed into the main body of the blue genes . after each resampling , we identify a cluster that contains more than 80% of all the blue genes , but as few non - blue genes as possible ( measured as a percentage of all genes in this cluster ) . this cluster is considered to be the cluster that corresponds to the original one in figure 8a . in figure 8b the root node of the ' corresponding cluster ' is marked with small horizontal line intersecting the vertical line ( representing the range of the cluster ) on the right of the clustering picture . then , for each of all the 225 genes , if it belongs to this ' corresponding cluster ' , we increase its ' in - cluster ' count by 1 . after resampling 30 times , the in - cluster counts are indicated in gray - scale on the left side of the original clustering ( figure 8c ) , with black representing 30 and white representing zero . a high ' in - cluster ' count indicates a gene ' remains ' in the original cluster in most of the resampled clustering trees . gene clustering ( a ) 225 filtered genes are clustered based on their expression profiles across 20 samples . each gene 's expression values are standardized to have mean 0 and sd 1 across 20 samples . although the original ' blue ' genes are scattered to various places , we can still determine where the original cluster is , using the criteria described in the text . ( c ) after resampling 30 times , the reliability of the genes belonging to the original cluster is indicated by the vertical gray - scale bar on the left of the blue - red picture . we can see from figure 8c that most genes in the original cluster are reliable members , whereas a few genes at the bottom of the cluster are not ( in fact they are merged into the original cluster last ) . interestingly , some genes originally not in the original cluster group with the ' corresponding clusters ' during resampling many times and have gray ' in - cluster ' marks . these genes may be related to the original cluster in some way . in summary , this method can help us to distinguish reliable and unreliable gene members of a cluster , as well as draw our attention to related genes originally clustered somewhere else because of the accidental nature of hierarchical clustering . we have developed a software package dna - chip analyzer ( dchip ) to perform invariant - set normalization ( see below ) , calculation of mbei , computation of confidence intervals of fold changes , and hierarchical clustering with resampling . our experience is that more than 10 arrays are appropriate for model training , outlier detection and mbei calculation . researchers with fewer than 10 arrays may seek arrays of the same chip type and hybridizing to similar tissue samples , and combine them in a single dchip analysis session . we are exploring model - based meta - analysis of many arrays of the same chip type but hybridizing to a heterogeneous set of tissues samples , and will present such analysis in future work . as array images usually have different overall image brightness ( figure 9a ) , especially when they are generated at different times and places , proper normalization is required before comparing the expression levels of genes between arrays . model - based expression computation requires normalized probe - level data ( from affymetrix 's dat or cel files ) . for a group of arrays , we normalize all arrays ( except the baseline array ) to a common baseline array having the median overall brightness ( as measured by the median cel intensity in an array ) . ( a ) the cel intensities ( see text ) of a pair of replicate arrays ( array 11 and 12 in array set 5 ) are plotted against each other . the baseline array 11 ( shown on the y - axis ) is not as bright as array 12 ( shown on the x - axis ) . the smoothing spline ( green curve ) deviates from the diagonal line y = x ( blue curve ) , indicating the need for normalization . ( b ) the same plot as ( a ) with superimposed circles representing the invariant set , on the basis of which a piecewise linear normalization relationship is determined ( black dotted line , whose y - coordinate is the normalized value of array 12 ) . the normalization curve is close to the smoothing spline curve in ( a ) as the two arrays are replicated arrays and all probes should be invariant . ( c ) after normalization ( y - axis is the baseline array 11 , and x - axis the normalized value of array 12 ) , the scatterplot centers around the diagonal line and the array 12 is adjusted to have the similar overall brightness as array 11 . ( d ) the q - q plot of probe intensities of array 11 and normalized array 12 shows the probes in the two sets have almost the same distribution . a normalization relation can be understood as a curve in the scatterplot of two arrays with the baseline array drawn on the y - axis and the array to be normalized on the x - axis . a straight line running through the origin is a multiplicative normalization method ( genechip 's scaling method ) , and a smoothing spline through the scatterplot can also be used ( figure 9a , also see ) . we should base the normalization only on probe values that belong to non - differentially expressed genes , but generally we do not know which genes are non - differentially expressed ( control or housekeeping genes may also be variable across arrays ) . nevertheless , we expect that a probe of a non - differentially expressed gene in two arrays to have similar intensity ranks ( ranks are calculated in two arrays separately ) . we use an iterative procedure to identify a set of probes ( called the invariant set ) , which presumably consists of points from non - differentially expressed genes ( figure 9b ) . specifically , we start with points of all pm probes ( about 140,000 for hu6800 array ) . if a point 's proportion rank difference ( prd , absolute rank difference in two arrays divided by n = 140,000 ) is small enough , it is kept for the new set . here the threshold of being small is prd < 0.003 when a points 's average intensity ranks in the two arrays is small and prd < 0.007 when it is large , accounting for fewer points at high - intensity range ; and the threshold is interpolated in between . we chose these parameters empirically to make the selected points in the invariant set thin enough to naturally determine a normalization relation . in this way we may obtain a new set of 10,000 points , and the same procedure is applied to the new set iteratively , until the number of points in the new set does not decrease anymore . a piecewise linear running median line figure 10 shows another pair of arrays where the normalization relationship is non - linear . similar plots as in figure 9 for arrays hybridized to two different samples ( array 24 and 36 of array set 5 ) . ( a ) cel intensities ; ( b ) same plot as in ( a ) with superimposed circles representing the invariant set ; ( c ) after renormalization ; ( d ) q - q plot of normalized probe intensities . note that the smoothing spline in ( a ) is affected by several points at the lower - right corner , which might belong to differentially expressed genes . the invariant set , on the other hand , does not include these points when determining the normalization curve , leading to a different normalization relationship at the high end . we have developed a software package dna - chip analyzer ( dchip ) to perform invariant - set normalization ( see below ) , calculation of mbei , computation of confidence intervals of fold changes , and hierarchical clustering with resampling . our experience is that more than 10 arrays are appropriate for model training , outlier detection and mbei calculation . researchers with fewer than 10 arrays may seek arrays of the same chip type and hybridizing to similar tissue samples , and combine them in a single dchip analysis session . we are exploring model - based meta - analysis of many arrays of the same chip type but hybridizing to a heterogeneous set of tissues samples , and will present such analysis in future work . as array images usually have different overall image brightness ( figure 9a ) , especially when they are generated at different times and places , proper normalization is required before comparing the expression levels of genes between arrays . model - based expression computation requires normalized probe - level data ( from affymetrix 's dat or cel files ) . for a group of arrays , we normalize all arrays ( except the baseline array ) to a common baseline array having the median overall brightness ( as measured by the median cel intensity in an array ) . ( a ) the cel intensities ( see text ) of a pair of replicate arrays ( array 11 and 12 in array set 5 ) are plotted against each other . the baseline array 11 ( shown on the y - axis ) is not as bright as array 12 ( shown on the x - axis ) . the smoothing spline ( green curve ) deviates from the diagonal line y = x ( blue curve ) , indicating the need for normalization . ( b ) the same plot as ( a ) with superimposed circles representing the invariant set , on the basis of which a piecewise linear normalization relationship is determined ( black dotted line , whose y - coordinate is the normalized value of array 12 ) . the normalization curve is close to the smoothing spline curve in ( a ) as the two arrays are replicated arrays and all probes should be invariant . ( c ) after normalization ( y - axis is the baseline array 11 , and x - axis the normalized value of array 12 ) , the scatterplot centers around the diagonal line and the array 12 is adjusted to have the similar overall brightness as array 11 . ( d ) the q - q plot of probe intensities of array 11 and normalized array 12 shows the probes in the two sets have almost the same distribution . a normalization relation can be understood as a curve in the scatterplot of two arrays with the baseline array drawn on the y - axis and the array to be normalized on the x - axis . a straight line running through the origin is a multiplicative normalization method ( genechip 's scaling method ) , and a smoothing spline through the scatterplot can also be used ( figure 9a , also see ) . we should base the normalization only on probe values that belong to non - differentially expressed genes , but generally we do not know which genes are non - differentially expressed ( control or housekeeping genes may also be variable across arrays ) . nevertheless , we expect that a probe of a non - differentially expressed gene in two arrays to have similar intensity ranks ( ranks are calculated in two arrays separately ) . we use an iterative procedure to identify a set of probes ( called the invariant set ) , which presumably consists of points from non - differentially expressed genes ( figure 9b ) . specifically , we start with points of all pm probes ( about 140,000 for hu6800 array ) . if a point 's proportion rank difference ( prd , absolute rank difference in two arrays divided by n = 140,000 ) is small enough , it is kept for the new set . here the threshold of being small is prd < 0.003 when a points 's average intensity ranks in the two arrays is small and prd < 0.007 when it is large , accounting for fewer points at high - intensity range ; and the threshold is interpolated in between . we chose these parameters empirically to make the selected points in the invariant set thin enough to naturally determine a normalization relation . in this way we may obtain a new set of 10,000 points , and the same procedure is applied to the new set iteratively , until the number of points in the new set does not decrease anymore . a piecewise linear running median line figure 10 shows another pair of arrays where the normalization relationship is non - linear . similar plots as in figure 9 for arrays hybridized to two different samples ( array 24 and 36 of array set 5 ) . ( a ) cel intensities ; ( b ) same plot as in ( a ) with superimposed circles representing the invariant set ; ( c ) after renormalization ; ( d ) q - q plot of normalized probe intensities . note that the smoothing spline in ( a ) is affected by several points at the lower - right corner , which might belong to differentially expressed genes . the invariant set , on the other hand , does not include these points when determining the normalization curve , leading to a different normalization relationship at the high end . we thank sven de vos , dan tang , nik brown , stan nelson , jae k. lee , yaron hakak , john walker and arindam bhattacharjee for providing data , and the editor and referees who provided valuable suggestions . this work is supported in part by nih grant 1 ro1 hg02341 - 01 and nsf grant dbi-9904701 .