Statistical outlier detection for gene expression microarray data

[0002] 1. Technical Field

[0003] The present invention is generally directed to the field of processing genomic data. More specifically, the invention relates to a system and method for performing statistical outlier detection for gene expression microarray data.

[0004] 2. Description of the Related Art

[0005] In genomics research, gene expression arrays are a breakthrough technology enabling the measurement of tens of thousands genes' transcription simultaneously. Because the numerical data associated with expression arrays usually arises from image processing, data quality is an important issue.

[0006] Two recent scientific articles, Schadt et al. (2000) and Li and Wong (2001), discuss this data quality issue for one of the most popular expression array platforms, the Affymetrix GeneChip™. For example, they point out that outlier problems may arise due to particle contaminations (see, FIG. 1 in Schadt et al. (2000)) or scratch contaminations (see FIG. 5 in Li and Wong (2001)). They indicate that improper statistical handling of aberrant or outlying data points can mislead analysis results.

[0007] Li and Wong propose an outlier detection method based on a multiplicative statistical model. While this approach is useful, it is limited to Affymetrix data and lacks the flexibility to accommodate more complex experimental designs. The multiplicative model used by the Li and Wong is as follows:

[0008] Y_ij=θ_iΦ_j+ε_ij, Σ_jΦ_j²=J, ε_ij˜N(0, σ²).  (1)

[0009] Y_ijis the intensity measurement of the j^thprobe in the i^tharray. θ_iis the i_thfixed array effect, Φ_jis the j^thfixed probe effect, and J is the number of probes. The ε_ij′s are assumed to be independent identically distributed normal random variables with mean 0 and variance σ². With the assumption of knowing Φs or θs, the following conditional means and standard errors can be derived and used in the Li and Wong method. $\begin{array}{c}{\tilde{\theta }}_i=\frac{{\sum }_j{Y}_{\mathrm{ij}}{\Phi }_j}{{\sum }_j{\Phi }_j^2},{\Phi }_j=\frac{{\sum }_i{Y}_{\mathrm{ij}}{\theta }_i}{{\sum }_i{\theta }_i},\\ \mathrm{StdErr}\left({\tilde{\theta }}_i\right)=\sqrt{\frac{{\sum }_j{\left(Y_{\mathrm{ij}}-{\hat{Y}}_{\mathrm{ij}}\right)}^2}{J\left(J-1\right)}},\\ \mathrm{StdErr}( {\tilde{\Phi }}_i)=\sqrt{\frac{{\sum }_i{\left(Y_{\mathrm{ij}}-{\hat{Y}}_{\mathrm{ij}}\right)}^2}{K\left(K-1\right)}},\\ K={\sum }_i{\tilde{\theta }}_i^2.\end{array}$

[0010] The following is a description of the Li and Wong outlier detection approach:

[0011] 1. Check array outliers—Fit the model (1) and calculate the conditional standard errors for all θ_i′s. Designate array as array outlier if either of the following criteria are met:

[0012] i. Associated θ has standard error larger than three times the median standard error of all θ_i′s.

[0013] ii. Associated θ has dominating magnitude with square value larger than 0.8 times the sum of squares of all θs.

[0014] Select out those array outliers and go to step 2.

[0015] 2. Check probe outliers—Fit the model (1) and calculate the conditional standard error for all Φ_j′s. Designate probe as probe outlier if either of the following criteria are met:

[0016] i. Associated Φ has standard error larger than three times the median standard error of all Φ_j′s.

[0017] ii. Associated Φ has dominating magnitude with square value larger than 0.8 times the sum of squares of all θ_j′s.

[0018] Select out those probe outliers and go to step 3.

[0019] 3. Iterate steps 1 and 2 until no further array or probe outliers selected.

[0020] In accordance with the disclosure below, a computer-implemented method and system are provided for detecting outliers in microarray data. A mixed linear statistical model is used to generate predictions based upon the received microarray data. Residuals are generated by subtracting model-based predictions from the original microarray sample data. Statistical tests are performed for residuals by adding covariates to the mixed model and testing their significance. Data from the microarrays are designated as outliers based upon the tested significance.

[0021] The present invention satisfies the general needs noted above and provides many advantages, as will become apparent from the following description when read in conjunction with the accompanying drawings, wherein:

[0022] FIG. 1 is a block diagram depicting the environment of the gene sample statistical modeling system;

[0023] FIGS. 2 and 3 are block diagrams depicting different components of the gene sample statistical modeling system;

[0024] FIG. 4 is a block diagram depicting components of the gene sample statistical modeling system used with additional statistical methods;

[0025] FIGS. 5A and 5B are flowcharts showing exemplary method steps and system elements for processing a GeneChip oligonucleotide microarray;

[0026] FIG. 6 is a block diagram depicting an example of a GeneChip probe design setting and experimental design layout;

[0027] FIGS. 7A and 7B are charts depicting data resulting from processing the GeneChip probe of FIG. 6;

[0028] FIG. 8 is a chart depicting mean values associated with the processing of the GeneChip probe of FIG. 6;

[0029] FIGS. 9A-9C show x-y graphs depicting expression profiles after using the statistical modeling system;

[0030] FIGS. 10A-10C show x-y graphs depicting expression profiles after using a prior art gene sample processing system;

[0031] FIG. 11 is a block diagram depicting use of a number of different arrays with the statistical modeling system;

[0032] FIG. 12 is a flowchart showing exemplary method steps and system elements for processing a cDNA microarray; and

[0033] FIG. 13 is a block diagram depicting use of the statistical modeling system with a number of applications.

[0034] FIG. 1 depicts a statistical modeling system 40 for use within a gene sample analysis environment. Experimenters may wish to analyze gene samples for many reasons, such as to better understand gene hybridization under different conditions or treatments. A preliminary step in the gene analysis process is to prepare the gene samples so that they may properly undergo hybridization. Preparation 30 may include attaching gene fragments (spots) onto glass slides so as to form miniature arrays.

[0035] The prepared gene samples are processed at 32 by one of many techniques. One technique involves hybridizing the gene samples, and then processing the gene samples to obtain image intensity data, such as by reading the intensity of each spot with a fluorescent detector (or other such device such as a charge-couple device). The image data from the processed samples 34 leads to gene-specific numerical intensities representing relative expression levels, and these in turn form the data set input 36 to computational analysis designed to assess relationships across biological samples.

[0036] Due to the reasons mentioned above, the data set 36 may contain outliers 38, and therefore it is desirable to have the outliers 38 eliminated before proceeding with the computational analysis. The statistical modeling system 40 is designed to identify and eliminate outliers 38 from the data set 36 and to generate a statistical model 60 that has better predictive capability due to the outliers being eliminated 62. FIG. 2 shows additional detail of the statistical modeling system 40 performing outlier detection and elimination.

[0037] With reference to FIG. 2, the statistical modeling system 40 utilizes one or more components 70 to identify and eliminate outliers 38 from the input data set 36. One component may involve the use of a mixed linear model 62. As used herein, the term “mixed” is broadly defined as containing both fixed and random factor effects. A factor effect is fixed if all the possible levels about which inference is to be made are represented in the study. A factor is random if the levels used in the study represent only a random sample of a larger set of potential levels. The statistical modeling system 40 may also perform rigid hypothesis testing 72 on putative outliers 38. Hypothesis testing occurs after an appropriate mixed model is fit to the data set and standardized residuals are calculated and ranked according their absolute magnitude. A statistical test for the residuals can be constructed by adding additional covariates to the model, each of which is an indicator variable for the observation in question. An indicator variable contains a 1 for the particular observation and a 0 for all other observations. The method then refits the model and tests the statistical significance of the covariates to decide whether the indicated observations are indeed outliers. If they are outliers, the statistical modeling system 40 eliminates the associated observations and searches for new outliers. If they are not outliers, the outlier analysis concludes.

[0038] The hypothesis testing operation 72 can be tuned 74 to select one or more outliers at a time and to test them for statistical validity. Tuning parameters are also available for determining the size of the initial group of potential outliers and for the selectivity of the statistical test.

[0039] It should be noted that the statistical modeling system 40 may use one or more of the components to identify and eliminate outliers 38 from the data set 36. For example, the residuals may be analyzed in a way other than through hypothesis testing 72, such as by directly eliminating the data points whose residuals have absolute magnitudes larger than a predetermined threshold such as three times the estimated standard deviation.

[0040] The components 70 may include additional capabilities, such as shown in FIG. 3 wherein the mixed model 70 includes effects with multiple degrees of freedom for a particular effect 80. Instead of previous approaches that are limited to a single degree of freedom (e.g., θ in model (1)), the mixed model 72 may have compound effects. For example, effects such as array, cell line, treatment, and cell line-treatment interaction may have more than two degrees of freedom. This allows greater robustness in model prediction of the data set 36 as well as an extension to accommodate different types of experimental designs (e.g., a circular design as described by Kerr and Churchill (2000), or an incomplete block and split plot design, or varying the level of replications, such as by spotting genes multiple times per array or using mRNA samples on multiple arrays, etc.).

[0041] Still further to illustrate the wide scope of the statistical modeling system 40, the statistical modeling system 40 may utilize one or more of the components 70 in combination with additional statistical approaches 90. As shown in FIG. 4, the statistical modeling system 40 may use a mixed model 70 approach (with subsequent residual analysis) in order to perform an initial outlier determination for the data set 36. The statistical modeling system 40 may then use another statistical approach 90, such as the above described Li and Wong approach (in whole or in part) to perform additional outlier analysis. It should be understood that the order of when to use other statistical approaches 90 (if at all) with respect to the components 70 may change based upon the application at hand. For example, another statistical approach 90 may first be used for outlier analysis, and then the components 70 of the statistical modeling system 40 may then provide additional outlier analysis

[0042] FIGS. 5A and 5B set forth a flowchart diagram 110 showing exemplary steps and system elements of the statistical modeling system in detecting outliers associated with a GeneChip oligonucleotide microarray. The method begins with the input data 112 being fit to an appropriate mixed linear model at process block 114. From here, standardized residuals are calculated in process block 116. The standardized residuals are then compared to a predetermined cutoff (c) in decision block 118. Because the standardized residuals are typically approximately normally distributed, the cutoff here can be determined accordingly. For example, setting the cutoff at 1.96 will examine 5% of data that exhibit extreme standardized residuals. If no standardized residual is greater than the predetermined cutoff, then processing continues at process block 142 where no outliers are designated, the results from the last time fitting template model are saved for further analysis at process block 132 and the method ends at termination block 134.

[0043] If a standardized residual is greater than the predetermined cutoff, then control passes to process block 120, where a number of outliers (n) are selected to be considered together and the model is refit with additional covariates corresponding to the indicator variables for the n observations with the largest standardized residuals in absolute value. The significances of the n covariates are tested to determine whether any calculated p-values are less than a predetermined significant level (a) at decision block 122. The model fitting and p-values can be calculated using a restricted maximum likelihood approach (see for a general discussion: Searle et al. (1993) Variance Components. Wiley, N.Y.). If the p-values are less than a, then control passes to process block 136 where observations with significant covariates are designated to be outliers and set to missing. Control passes back to process block 114, and the template model is refit with the modified data.

[0044] If the p-values are greater than a, control passes to process block 124 where the standard deviation of the standardized residual for each probe is calculated (PSD) across arrays and the standard deviation of PSD (SDPSD) are calculated across probes. Decision block 126 then tests whether any PSD exceeds three times the SDPSD from the average of the PSDs. If a PSD does exceed that amount, then process block 138 sets the probes to missing and control is returned to process block 114. If no PSDs exceed three times the SDPSD from the average of the PSDs, then control passes to process block 128.

[0045] In process block 128, a standard deviation of standardized residuals for each array (ASD) is calculated. Then, the standard deviation of ASD (SDASD) across arrays is calculated. Decision block 130 tests whether any ASD exceeds three times the SDASD from the average of the ASD. If any ASD exceeds three times the SDASD from the average of the ASD, then process block 140 sets the array to missing and returns control to process block 114. If no ASD exceeds three times the SDASD from the average of the ASD, then control passes to process block 132 where the results from the previous fitted template model are saved for further analysis and the method ends at termination block 134.

[0046] FIG. 6 illustrates an exemplary GeneChip probe design setting and experimental design layout. In this example, there are eight GeneChips (150, 152, 154, 156, 158, 160, 162, and 164). Two treatments 170 and 172 were applied on two cell lines with two replicates. P1-P20 are the twenty probe intensity measurements per chip. The experiment studied the gene expression of each probe under an irradiated treatment condition and an unradiated condition. It should be understood that many different array configurations and experimental design layouts may be used and the instant configuration of FIG. 6 is only an example.

[0047] The data for this example constitute artificial expression values for one gene based on the experimental design of the ionizing radiation response data used in Tusher et al. (2001). The data set is listed in FIGS. 7A and 7B. With reference to FIGS. 7A and 7B, chip identification information is provided in column 200; line identification information is provided in column 202; treatment identification information is provided in column 204; replicate identification information is provided in column 206; probe identification information is provided in column 208; and column 210 indicates the intensity values obtained from processing the chips. As an illustration, probe #1 in line #1 of chip #1 underwent an irradiation treatment and provided an intensity value of 744.8. Of note are the two entries for probes #11 and #12 of chip #3. These are true outliers in the data set which are to be correctly identified and eliminated by the statistical modeling system while not eliminating non-outlier points from the data set.

[0048] The mixed model used in this example is as follows: $\begin{array}{cc}{\log }_2\left(Y_{\mathrm{ijkl}}\right)-M_{\mathrm{ijl}}=L_i+T_j+{\mathrm{LT}}_{\mathrm{ij}}+P_k+{\mathrm{LP}}_{\mathrm{ik}}+{\mathrm{TP}}_{\mathrm{jk}}+A_{l\left(\mathrm{ij}\right)}+{ϵ}_{\mathrm{ijkl}}.,A_{l\left(\mathrm{ij}\right)}+{ϵ}_{\mathrm{ijkl}}\sim N\left(0,{\sigma }_a^2+{\sigma }^2\right),\mathrm{Cov}\left(A_{l\left(\mathrm{ij}\right)}+{ϵ}_{\mathrm{ijkl}},A_{l^{\prime }\left(i^{\prime }{j}^{\prime }\right)}+{ϵ}_{i^{\prime }{j}^{\prime }{k}^{\prime }{l}^{\prime }}\right)=\{\begin{array}{cc}{\sigma }_a^2+{\sigma }^2& \mathrm{if}\left(i,j,k,l\right)=\left(i^{\prime },j^{\prime },k^{\prime },l^{\prime }\right),\\ {\sigma }_a^2& \mathrm{if}\left(i,j,l\right)=\left(i^{\prime },j^{\prime },l^{\prime }\right)\mathrm{but}k\ne k^{\prime },\\ 0& \mathrm{otherwise}.\end{array}& \left(2\right)\end{array}$

[0049] Y_ijklis the expression measurement of the i^thcell line applying the j^thtreatment at the k^thprobe in the l^threplicate. M_ijlis the mean value of logged (base 2) intensity of the chip associated with the i^thcell line, the j^thtreatment, and the l^threplicate across all genes. For this example, the eight intensity mean values are listed in FIG. 8 and used for data normalization prior to modeling. The symbols L, T, LT, P, LP, TP and A represent cell line, treatment, cell line-treatment interaction, probe, cell line-probe interaction, treatment-probe interaction, and array effects, respectively. The A_l(ij)′s are assumed to be independent and identically distributed normal random variables with mean 0 and variance σ_a². The ε_ijkl′s are assumed to be independent identically distributed normal random variables with mean 0 and variance σ², and are independent of the A_l(ij)′s. The remaining terms in the model are assumed to be fixed effects. The use of both fixed and random effects in the model allow the model to be a mixed model.

[0050] In this example, there are three tuning parameters, c, n, and a. The parameter c is used to ensure that the prospective outliers are far enough from the prediction of the model. We set c=1.645 in the example, which according to a normality assumption will result in checking the 10% most extreme standardized residuals. The parameter n is the number of outliers to consider at one time; typically n=1. The parameter a is used as the significance level for testing indicator-variable covariates. For this example, we set a=0.0026 for subsequent comparisons to the Li and Wong method because their “three standard error rule” roughly controls at the 0.0026 significance level if the standardized residuals are independent and identically distributed normal random variables. It should be understood that different parameter values may be used depending upon the application at hand.

[0051] The plots 250 of FIGS. 9A-9C show the expression profiles of the gene in all eight arrays using the statistical modeling system (which is labeled in this example as “FOD” (forward outlier detection method)). The three curves in each plot 250 represent the original measurements (Y), predictions from the mixed model (Mix) after first time fitting, and predictions after the method has eliminated outliers (FOD). For example in plot 252, curve 260 depicts the original measurements (Y) for array “1un1” (i.e., the array labeled as replicate #1, line #1, with the unradiated treatment); curve 262 depicts the predictions from the mixed model (MIX); and curve 264 depicts the predictions after the method has eliminated outliers (FOD). After comparison of the original measurement curves 260 for these eight arrays, the measurements 270 and 272 of probes 11 and 12 in array “1un1” in plot 252 are determined by the FOD method to be outliers whereas the measurements 280 and 282 of plot 254 are determined not to be outliers.

[0052] Note that if the original measurement curves 260 and mix curves 262 for plots 252 and 254 are compared, then both probe 11 at 270 in plot 252 and probe 11 at 280 in plot 254 will be selected as outliers if selecting multiple outliers at each time. In this case, the FOD method only selects probes 11 and 12 in array 1un1 as single-point outliers.

[0053] The plots 300 of FIGS. 10A-10C show the expression profiles of the gene in the eight arrays described in FIGS. 9A-9C. The expression profiles of FIGS. 10A-10C were generated using model (1) and the outlier detection approach disclosed in Li and Wong (2001). The three curves in each plot 300 represent the original measurements (Y), predictions from model (1) (MuP1) after first time fitting, and predictions after applying the LW outlier detection approach (LW). For example in plot 302, curve 310 represents the original measurements (Y); curve 312 represents predictions from model (1) (MuP1) after first time fitting; and curve 314 represents predictions (LW) after applying the LW outlier detection approach.

[0054] In this case, the LW approach selects no array outlier but “probe 11” 320 of all eight arrays (318, 320, 324, 326, 328, 330, 332, 334) as probe outliers in the first iteration. After eliminating “probe 11” 320, the “probe 12” 322 makes the within array standard error of array “1un1” 304 about six times larger than the median within array standard error. Therefore, the LW approach further selects the entire array “1un1” 304 as array outliers in the second iteration. In total, the LW approach selects twenty-seven observations as outliers. However, note that if “probe 11” 320 and “probe 12” 324 in array “1un1” 304 are removed from the data set, then the LW approach does not select any outliers. Thus, under model (1) and the LW approach, “probe 11” 320 and “probe 12” 324 in array “1un1” 304 are two extremely influential observations which cause the incorrect classification of the other twenty-five observations as outliers by this approach. This may be due to the LW approach not examining single-point outliers first, but rather applying the “3 standard error rule” and then conservatively selecting array and probe outliers (which actually result from few of the single-point outliers as described in the example). Table 1 highlights additional differences between FOD and LW.

[0055] The preferred embodiment described with reference to the drawing figures is presented only to demonstrate examples of the invention. Additional and/or alternative embodiments of the invention should be apparent to one of ordinary skill in the art upon reading this disclosure as the present invention is applicable in many contexts. As an example of the wide scope of the invention, FIG. 11 depicts the statistical modeling system 40 detecting outliers in data arising from any kind of gene expression microarray platform (such as two-color spotted arrays 400, multi-color spotted arrays 402, nylon filter arrays 404, etc.).

[0056] For example, the statistical modeling system and method may process data from cDNA microarrays. Processing for of this type of microarray is shown by the flowchart 450 of FIG. 12. With reference to FIG. 12, the method begins with the input data 452 being fit to an appropriate mixed linear model at process block 454. From here, standardized residuals are calculated in process block 456. The standardized residuals are then compared to a predetermined cutoff (c) in decision block 458. Because the standardized residuals are distributed by standard normal, the cutoff here may be determined according to this normality. For example, setting the cutoff at 1.96 will examine 5% of data that exhibit extreme standardized residuals. If no standardized residual is greater than the predetermined cutoff, then processing continues at process block 470 where no outliers are designated, the results from the previous fitted template model are saved for further analysis at process block 464 and the method ends at termination block 466.

[0057] If a standardized residual is greater than the predetermined cutoff, then control passes to process block 460, where a number of outliers (n) are selected to be considered together and the model is refit with additional covariates corresponding to the indicator variables for the n observations with the largest standardized residuals in absolute value. The significance of the n covariates are tested to determine whether any calculated p-values are less than a predetermined p-value (a) at decision block 462. If the p-values are less than a, then control passes to process block 468 where observations with significant covariates are set to missing. Control passes back to process block 454 and a new model is fit to the remaining covariates. If the p-values are greater than a, then the results from the previous fitted template model are saved for further analysis at process block 464 and the method ends at termination block 466.

[0058] As yet another example of the wide scope of the system and method, FIG. 13 shows use of the statistical modeling system with a number of different applications 500. For example, the data model 60 (with outliers eliminated) may be used in any subsequent analysis, such in cluster analysis. More specifically, experimenters may use the data model 60 as a precursor to clustering to ensure the inputs are statistically meaningful, or it can be used after clustering to explore and validate implied associations. It should be understood that not only may the generated data models be used by different applications, but also the original data set (with outliers eliminated) may be used.