A robust statistical procedure to discover expression biomarkers using microarray genomic expression data

A robust statistical procedure to discover expression biomarkers using microarray genomic expression data first-author corresp ‡Corresponding Author ‡Corresponding Author †E-mail:jzhu@zju.edu.cn †E-mail:jzhu@zju.edu.cn jzhu@zju.edu.cn Institute of Bioinformatics, Zhejiang University, Hangzhou 310029, China Institute of Bioinformatics, Zhejiang University, Hangzhou 310029, China Microarray has become increasingly popular biotechnology in biological and medical researches, and has been widely applied in classification of treatment subtypes using expression patterns of biomarkers. We developed a statistical procedure to identify expression biomarkers for treatment subtype classification by constructing an F -statistic based on Henderson method III. Monte Carlo simulations were conducted to examine the robustness and efficiency of the proposed method. Simulation results showed that our method could provide satisfying power of identifying differentially expressed genes (DEGs) with false discovery rate (FDR) lower than the given type I error rate. In addition, we analyzed a leukemia dataset collected from 38 leukemia patients with 27 samples diagnosed as acute lymphoblastic leukemia (ALL) and 11 samples as acute myeloid leukemia (AML). We compared our results with those from the methods of significance analysis of microarray (SAM) and microarray analysis of variance (MAANOVA). Among these three methods, only expression biomarkers identified by our method can precisely identify the three human acute leukemia subtypes. Microarray Biomarker Henderson method III Gene expression pattern Mixed linear model INTRODUCTION 1 Microarray technique is a powerful laboratory tool for simultaneously monitoring genome-wide expression in different conditions. One of the most important applications of microarray technique is the classification of tumor subtypes or different disease states to facilitate clinical researchers in diagnostic, therapeutic or prognostic decisions for patients (Golub et al., 1999 1999 ; Alizadeh et al., 2000 2000 ; Spindler, 2006 2006 ). The generally used approaches, such as cluster analysis and supervised grouping, only focus on the similarity of the data structure, and fail to guarantee that the used class predictors are biologically associated with class distinction. Therefore, prior to cluster analysis, a central need is to explore whether there are some expression biomarkers strongly correlated with specific classes. Some statistical methods such as significance analysis of microarray (SAM) (Tusher et al., 2001 2001 ), fixed ANOVA method (Kerr et al., 2000 2000 ; Kerr and Churchill, 2001 2001 ) and mixed linear model method (Wolfinger et al., 2001 2001 ; Jin et al., 2001 2001 ; Lu et al., 2005 2005 ) have been proposed to serve this goal. Based on Henderson method III, we developed a statistical method under the mixed linear model framework (Zhu, 2000 2000 ) to objectively identify expression biomarkers for treatment classification. Monte Carlo simulations were conducted to validate the robustness and efficiency of the present method, and a real dataset of leukemia was analyzed to assess the utility of the method. METHOD FRAMEWORK 1 As in most analysis methods of microarray, we first use the normalization procedure to minimize the global systematical variations involved in the experiment from the original fluorescence measurements. The normalization model can be written as y ijklm = μ + A i + D j + P k + T l + γ ijklm , where μ represents the mean expression level over all genes, fixed effect; A i is the array effect, random effect, $$A_i \sim(0,\sigma _A^2 )$$ ; D j is the dye effect, random effect, $$D_j \sim(0,\sigma _D^2 )$$ ; P k is the pin effect, random effect, $$P_k \sim(0,\sigma _P^2 )$$ ; γ ijklm is residual error, $$\gamma _{ijklm} \sim(0,\sigma _\gamma ^2 )$$ . γ ijklm is obtained by subtracting the fitted values of the effects in model (1) (1) from base 2 logarithm of background-corrected measurements ( y ijklm ) using least square estimation (LSE) method, and will subsequently be used as the inputs for the gene-specific models γ ijglm = μ g + A gi + D gj + T gl + ε ijglm , where μ g represents the mean expression level of gene g ; A gi is gene specific array effect, $$A_{gi} \sim(0,\sigma _{A_g }^2 )$$ ; D gj is gene specific dye effect, $$D_{gj} \sim(0,\sigma _{D_g }^2 )$$ ; ε ijglm is gene-dependent residual error, $$\varepsilon _{ijglm} \sim(0,\sigma _{\varepsilon _g }^2 )$$ , which is different from γ ijklm in model (1) (1) . Under the null hypothesis H 0 : T g 1 = T g 2 =…=0, Henderson method III is employed to construct the F -statistic to test the significance of treatment effects (Searle, 1971 1971 ). Since analysis of microarray data involves multiple statistical tests, we use false discovery rate (FDR) (Benjamini and Hochberg, 1995 1995 ) to control the experimental-wise type I error. The identified differentially expressed genes (DEGs) are ranked by their statistic scores which can provide more information and choice for biologists. Finally, the potentially DEGs detected by model (2) (2) are fitted in the following full model to estimate the variance components and effects interested. y ijkglm = μ + G g + A i + D j + P k + T l + GA gi + GD gj + GT gl + ε ijkglm , where μ and G g are fixed effects, and A i , D j , P k , GA gi , GD gj , ε ijkglm are random effects normally distributed with zero means and variance components σ A 2 , σ D 2 , σ P 2 , σ GA 2 , σ GD 2 , σ ε 2 , respectively. The terms T l , T gl and GT gl in models can be regarded as fixed effects or random effects according to the experimental intent. Since different pins have different characteristics and surface properties with different amounts of target cDNA, we include the P k effect in our model. These models are extensible to more complex situations such as N -dyes, multiple factors decomposed from the treatment effect and other variations like fluctuations due to mRNA extraction, cDNA synthesis. Markov Chain Monte Carlo (MCMC) method (Wang et al., 1994 1994 ) is used to estimate the variance components of random effects in the model, to estimate fixed effects, and to predict random effects as well. RESULTS 1 Simulation analysis 2 Different variation magnitudes were set according to the results from the previously analyzed real dataset available in Stanford microarray database. We assigned the residual variance as 1 with the proportion of gene by treatment interaction effect ( V GT ) variance to the residual variance ( V GT / V ε ) ranging from 1 to 10. Two assumptions were adopted: (1) All parameters in the model follow independent and identical normal distribution, denoted as NormAspt; (2) The observations from different array do not share the same variance, denoted as ArrayHetero. Simulation datasets were generated from different experimental designs, loop design with spots replicated within single array (denoted as LOOPREPIN), loop design with spots replicated between arrays (denoted as LOOPREPOUT), reference design with spots replicated within single array (denoted as REFREPIN), reference design with spots replicated between arrays (denoted as REFREPOUT). All the simulation datasets were run with 200 replicates, with powers of identifying DEGs with FDR control at 0.05 being shown in Fig. 1 1 . Simulation results revealed that our method was appealing for identifying DEGs validated with high power when V GT / V ε exceeded 2, especially in the case of LOOPREPIN. In addition, it was shown that our method could offer intriguing stability under different assumptions, which would be important for our method for analyzing microarray data when array variance heterogeneity is quite common in microarray data. However, in the case of REFREPOUT, the result was not perfect, due to the severe confounding between other variations like GA or GD and variance of GT effect in this experimental design. Therefore, it is strongly recommended that appropriate experimental design should be used, such as loop design with spots replicated within an array or multi-color microarray design (Woo et al., 2005 2005 ). Worked example 2 Previous study of leukemia (Golub et al., 1999 1999 ) monitored expression patterns from 38 leukemia patients (with clinically predefined of T-cell ALL, B-cell ALL and AML) to develop an expression-based molecular classification method for acute leukemia as an assistant tool of clinical diagnosis. Affymetrix Hu6800 GeneChips were used. Since the sample grouping of these datasets has been clinically verified, we could use it to validate the utilization of our method. And we also used the methods of SAM and MAANOVA to analyze this dataset for comparison. By a default in the configurations, MAANOVA could only identify 102 marker genes, less than those identified by the other two methods. Thus, we used the top ranked 102 genes from our method and SAM, as well as 102 significant genes from MAANOVA to classify the samples by hierarchical cluster using Pearson correlation distance with UPGMA-linkage criterion. In distinguishing two predefined classes of leukemia ALL and AML, our method yielded accurate classifications, while SAM confused to classify 3 ALL samples into the AML samples and MAANOVA incorrectly placed 8 ALL samples and 5 AML samples to each other. Expression biomarkers identified by our method were also sensitive for partitioning ALL samples into T-cell ALL and B-cell ALL subclasses with only two B-cell samples and one T-cell sample wrongly classified, while SAM and MAANOVA entirely failed to do it. Since putative cluster labels are available for leukemia data, external indices of adjusted rank index (Hubert and Arabie, 1985 1985 ), Jaccard index (Jain and Dubes, 1988 1988 ) and FM index (Fowlkes and Mallows, 1983 1983 ) were computed to evaluate the quality of cluster results. These external indices have the property that the higher the score, the better the cluster solution, with a score of 1.0 indicating a perfect solution. Compared cluster results showing dissimilar structure of expression data were due to the different biomarkers identified by our method, SAM, and MAANOVA, respectively, were summarized in Table 1 1 . Our method showed the highest scores in these three cluster validation measurements, indicating that the biomarkers discovered by our method were very close to these different classes of leukemia. Note: Two classes: ALL and AML subclasses; Three classes: T-cell ALL, B-cell ALL and AML 2 Method 3 Two classes 3 Three classes Adjusted rank Jaccard FM Adjusted rank Jaccard FM Our method 1.000 1.000 1.000 0.776 0.747 0.856 SAM 0.893 0.911 0.954 0.533 0.568 0.731 MAANOVA 0.012 0.383 0.554 0.076 0.252 0.403 DISCUSSION 1 The recognition of objective expression biomarkers plays a crucial role in correct classification of tumor subtypes which is valuable for assisting in clinical diagnosis. Many standard statistical methods have been used to address the issue, but none has yet obtained widespread acceptance because of the high rates of false discovery. In the present study, we implement a novel statistical approach in three interconnected steps, normalization (model (1) (1) ), gene-specific model fitting (model (2) (2) ) and multiple genes model fitting (model (3) (3) ). In the second step (model (2) (2) ), an F -statistic is constructed via Henderson method III to scale the expression change among different treatments of each gene. This strategy is quite efficient in terms of statistical power and computation. Simulation results in multiple configurations and the real leukemia data analysis showed that our method can improve the ability to correctly identify DEGs or expression biomarkers in expression profiles analysis. Meanwhile, technical and stochastic variations such as mRNA extraction, cDNA synthesis, labelling reactions and print or hybridization efficiency, are usually involved in microarray experiments, and inevitably lead to the noise in raw expression measurements and bias interpretation of class distinction. So it is recommended to do the third step—multi-gene model fitting (model (3) (3) ) using the MCMC method which can give unbiased prediction of GT interaction effects for cluster analysis or discriminant analysis. Besides, the estimates of various sources of variation can provide some feedback on the quality of the experiment to researchers, which is essential for improving the laboratory protocols for further experiments. For example, if variations of A and GA effects is large, it is essential to re-select appropriate array with inherently less variations, or use finer experimental design so as to construct appropriate statistical model to screen these variations. This is only a preliminary study, with the detailed research carried out in our following paper. References AA Alizadeh , MB Eisen , RE Davis , C Ma , IS Lossos , A Rosenwald , JC Boldrick , H Sabet , T Tran , X Yu , 等 . Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling . Nature , 2000 . 403 ( 6769 ): 503 - 511 . 10676951 DOI: 10.1038/35000501 http://doi.org/10.1038/35000501 . Distinct types of diffuse large B-cell lymphoma identified by gene expression profiling Y Benjamini , Y Hochberg . Controlling the false discovery rate: a practical and powerful approach to multiple testing . Journal of the Royal Statistical Society: Series B (Methodological) , 1995 . 57 ( 1 ): 289 - 300 . . Controlling the false discovery rate: a practical and powerful approach to multiple testing EB Fowlkes , CL Mallows . A method for comparing two hierarchical clusterings . J American Statistical Association , 1983 . 78 ( 383 ): 553 - 569 . DOI: 10.2307/2288117 http://doi.org/10.2307/2288117 . A method for comparing two hierarchical clusterings TR Golub , DK Slonim , P Tamayo , C Huard , M Gaasenbeek , JP Mesirov , H Coller , ML Loh , JR Downing , MA Caligiuri , 等 . Molecular classification of cancer: class discovery and class prediction by gene expression monitoring . Science , 1999 . 286 ( 5439 ): 531 - 537 . 10521349 DOI: 10.1126/science.286.5439.531 http://doi.org/10.1126/science.286.5439.531 . Molecular classification of cancer: class discovery and class prediction by gene expression monitoring L Hubert , P Arabie . Comparing partitions . Journal of Classification , 1985 . 2 ( 1 ): 193 - 218 . DOI: 10.1007/BF01908075 http://doi.org/10.1007/BF01908075 . Comparing partitions AK Jain , RC Dubes . Algorithms for Clustering Data , 1988 . : Englewood Cliffs, NJ Prentice Hall , . Algorithms for Clustering Data W Jin , RM Riley , RD Wolfinger , KP White , G Passador-Gurgel , G Gibson . The contributions of sex, genotype and age to transcriptional variance in drosophila melanogaster . Nature Genetics , 2001 . 29 ( 4 ): 389 - 395 . 11726925 DOI: 10.1038/ng766 http://doi.org/10.1038/ng766 . The contributions of sex, genotype and age to transcriptional variance in drosophila melanogaster MK Kerr , GA Churchill . Experimental design for gene expression microarrays . Biostatistics , 2001 . 2 ( 2 ): 183 - 201 . 12933549 DOI: 10.1093/biostatistics/2.2.183 http://doi.org/10.1093/biostatistics/2.2.183 . Experimental design for gene expression microarrays MK Kerr , M Martin , GA Churchill . Analysis of variance for gene expression microarray data . Journal of Computational Biology , 2000 . 7 ( 6 ): 819 - 837 . 11382364 DOI: 10.1089/10665270050514954 http://doi.org/10.1089/10665270050514954 . Analysis of variance for gene expression microarray data Y Lu , J Zhu , P Liu . A two-step strategy for detecting differential gene expression of cDNA microarray data . Current Genetics , 2005 . 47 ( 2 ): 121 - 131 . 15688252 DOI: 10.1007/s00294-004-0551-3 http://doi.org/10.1007/s00294-004-0551-3 . A two-step strategy for detecting differential gene expression of cDNA microarray data SR Searle . Linear Models , 1971 . : New York John Wiley & Sons , . Linear Models SR Spindler . Use of microarray biomarkers to identify longevity therapeutics . Aging Cell , 2006 . 5 ( 1 ): 39 - 50 . 16441842 DOI: 10.1111/j.1474-9726.2006.00194.x http://doi.org/10.1111/j.1474-9726.2006.00194.x . Use of microarray biomarkers to identify longevity therapeutics VG Tusher , R Tibshirani , G Chu . Significance analysis of microarrays applied to the ionizing radiation response . Proc Nat Acad Sci USA , 2001 . 98 ( 9 ): 5116 - 5121 . 11309499 DOI: 10.1073/pnas.091062498 http://doi.org/10.1073/pnas.091062498 . Significance analysis of microarrays applied to the ionizing radiation response CS Wang , JJ Rutledge , D Gianola . Bayesian analysis of mixed linear models via Gibbs sampling with an application to litter size in Iberian pigs . Genetics Selection Evolution , 1994 . 26 91 - 115 . . Bayesian analysis of mixed linear models via Gibbs sampling with an application to litter size in Iberian pigs RD Wolfinger , G Gibson , ED Wolfinger , L Bennett , H Hamadeh , P Bushel , C Afshari , RS Paules . Assessing gene significance from cDNA microarray expression data via mixed models . Journal of Computational Biology , 2001 . 8 ( 6 ): 625 - 638 . 11747616 DOI: 10.1089/106652701753307520 http://doi.org/10.1089/106652701753307520 . Assessing gene significance from cDNA microarray expression data via mixed models Y Woo , W Krueger , A Kaur , G Churchill . Experimental design for three-color and four-color gene expression microarrays . Bioinformatics , 2005 . 21 ( Suppl. 1 ): i459 - i467 . 15961491 DOI: 10.1093/bioinformatics/bti1031 http://doi.org/10.1093/bioinformatics/bti1031 . Experimental design for three-color and four-color gene expression microarrays J Zhu . Mixed linear model approaches for analyzing genetic models of complex quantitative traits. . Journal of Zhejiang University SCIENCE , 2000 . 1 ( 1 ): 78 - 90 . . Mixed linear model approaches for analyzing genetic models of complex quantitative traits. We thank Guobo Chen and Lide Han for providing helpful discussion and comments on this manuscript. 10.1631/jzus.2006.B0603 Project partly supported by the National Basic Research Program (973) of China (No. 2004CB117306) and the National Natural Science Foundation of China (No. 2002AA234031) 2006-04-04 2006-05-31 2006-07-14 2022-01-05 Copyright © 2006, Journal of Zhejiang University Science 2006 Journal of Zhejiang University-SCIENCE B (Biomedicine and Biotechnology) 08 07