Optimal Selection of Microarray Analysis Methods Using a Conceptual Clustering Algorithm

Optimal Selection of Microarray Analysis Methods Using a Conceptual Clustering Algorithm C. Rubio-Escudero1, R. Romero-Záliz1, O. Cordón1, O. Harari1, C. del Val1, and I. Zwir1,2 1 Department of Computer Science and Artificial Intelligence, C/Daniel Saucedo Aranda s/n, Granada 18071, Spain {crubio, rocio, ocordon, oharari, delval, zwir}@decsai.ugr.es 2 Howard Hughes Medical Institute, Washington University School of Medicine, St. Louis, MO Abstract. The rapid development of methods that select over/under expressed genes from microarray experiments have not yet matched the need for tools that identify informational profiles that differentiate between experimental conditions such as time, treatment and phenotype. Uncertainty arises when methods devoted to identify significantly expressed genes are evaluated: do all microarray analysis methods yield similar results from the same input dataset? do different microarray datasets require distinct analysis methods?. We performed a detailed evaluation of several microarray analysis methods, finding that none of these methods alone identifies all observable differential profiles, nor subsumes the results obtained by the other methods. Consequently, we propose a procedure that, given certain user-defined preferences, generates an optimal suite of statistical methods. These solutions are optimal in the sense that they constitute partial ordered subsets of all possible method-associations bounded by both, the most specific and the most sensitive available solution. 1 Introduction Advances in molecular biology and computational techniques permit the systematical study of molecular processes that underlie biological systems [1]. Particularly, microarray technology has revolutionized modern biomedical research by its capacity to monitor changes in RNA abundance for thousands of genes simultaneously [2]. To address the statistical challenge of analyzing these large data sets, new methods have emerged ([3], [4], [5], [6], [7] and many others). However, there is a dearth of computational methods to facilitate understanding of differential gene expression profiles (e.g., profiles that change over time and/or over treatments and/or over patient) and to decide which is the most reliable method to identify differences across profiles. We investigated the performance of several commonly used statistical methods, including T-Tests [4], Permutation Tests [5], Analysis of Variance [6] and Repeated Measures ANOVA [7], in identifying differential expression profiles that change over time, treatments and phenotype. We found that these methods do not identify all ob- servable distinct profiles. Moreover, none of them subsumes the results obtained by the other methods. In view of these results, we propose a conceptual clustering method [8], [9], [10], devoted to discover optimal associations of microarray analysis methods in an effort to identify differential gene expression profiles. 2 Methods We propose a conceptual clustering approach [8], [9], [10] devoted to identify optimal associations among microarray analysis methods in an effort to identify differential expression profiles (Fig. 1). This approach consists of six phases: (1) preprocessing of the dataset; (2) identification of differentially expressed genes by application of several statistical methods; (3) arrangement of a lattice structure containing all possible associations of the statistical methods applied; (4) association of differentially expressed genes into differential profiles by clustering genes that change their expression over time, patient and/or treatment; (5) evaluation of the performance of the method-associations based on their specificity and sensitivity in the identification of previously detected differential profiles, using multiobjective optimization techniques [11], [12]. We create a set of method association rules based on the learned mappings of differential profiles into method-associations, [13]; (6) finally, we are able to predict optimal method-associations to identify differential profiles in new microarray datasets by use of the method association rules. 2.1 Identification of Differentially Expressed Genes We perform the retrieval of differentially expressed genes from one experimental condition to the other/s by application of several statistical techniques [3], [14], harboring Student’s T-Test proposed in [4], including some of the variants the method poses to distinguish changes in the abundance of RNA occurring over both treatment and time; Permutation Test described in [5], also including a time approach; Analysis of Variance described in [6]; and Longitudinal Data approach by using Repeated Measures Analysis of Variance described in [7]. DIFFERENTIALLY EXPRESSED GENES MICROARRAY RAW DATA MICROARRAY PREPROCESSED DATA PREPROCESSING (1) (4) IDENTIFICATION OF DIFFERENTIAL PROFILES DIFFERENTIAL EXPRESSION PROFILES METHOD EVALUATION SCALING (2) IDENTIFICATION OF DIFF. EXPRESSED GENES (6) PREDICTION QUERY PROFILES + MICROARRAY DATA (OPTIONAL) METHOD SELECTION (5) CREATION OF METHOD ASS. RULES METHOD ASSOCIATION RULES (3) ASSOCIATION OF STATISTICAL METHODS LATTICE OF METHOD ASSOCIATIONS STATISTICAL METHODS CONCEPTUAL CLUSTERING PROFILE IDENTIF. PROFILE PRUNING NORMALIZING Fig. 1. Graphical representation of the methodology. The squared boxes represent the phases of the methodology, the round cornered boxes correspond to the input/output data at each step, and the ellipses the operations performed at each phase. 2.2 Detection of Method-Associations We arrange a lattice containing all potential associations of the statistical methods used to retrieve differentially expressed genes (Fig. 2). The methods are associated as: }...,...,,,,,{ 21312121 nn MMMMMMMMMMM ⊕⊕⊕⊕⊕= , (1) where ⊕ is a classical set operator (e.g., the union ( U) or the intersection (I)) applied to the sets of genes retrieved by each method, and 1 M corresponds to T-Test, 2 M to T-Test considering time, 3 M Permutation Test, 4 M Permutation Test con sidering time, 5 M ANOVA over treatment, 6 M ANOVA over time, 7 M ANOVA over treatment and time, 8 M RMANOVA over treatment, 9 M RMANOVA over time and 10 M RMANOVA over treatment and time. The lattice containing all potential method-associations, M, is structured from top (i.e., intersection of all methods) to bottom (i.e., union of all methods) [15]. Each node in the lattice ( MM i∈) is applied to the microarray dataset (D) retrieving the set of differentially expressed genes that are recognized by the method or methodassociations in such node D))Mi((. M 3 M 1  M 2  M 3   M n M 1 M 2 M 3 M 2 M 3 M 4 M n-2 M n-1 M n M 1 M 2 M 3 …M k M n-2 M n-1 M n M 2 M 3 M 4 M 1 M 2 M 3 M n-1 M n M 2 M 3 M 1 M 3 M 1 M 2 M n M 1 M 2 M 1 M 2 M 1 M 3 M 2 M 3 M n-1 M n M 3 …… …… … …… …… Fig. 2. Lattice structure containing all statistical methods potential associations 2.3 Identification of Differential Profiles The set of genes previously identified in Section 2.2 serves as a means to create differential expression profiles (i.e., sets of genes with coordinate changes in RNA abundance) between treatment T P, control C P and subject. The applied representation (Fig. 3) allows us to identify different pattern behavior among patients inside the same experimental group, since this information may be missed if patients in the same experimental group were not plotted individually. We clustered separately genes in treatment and control groups. Therefore, genes belonging to a cluster in treatment, T P, can fit in more than one cluster in control, C P, and vice versa. We apply the K-means clustering algorithm [16] and identify differential profiles denoted as )( CT PP , which are pairwise relationships between profiles, T P Treatment - 6 Control - 3 TREATMENT CONTROL HOUR 0 2 4 6 9 24 0 2 4 6 9 24 0 2 4 6 9 24 0 2 4 6 9 24 PATIENT 1 2 3 4 HOUR 0 2 4 6 9 24 0 2 4 6 9 24 0 2 4 6 9 24 0 2 4 6 9 24 PATIENT 5 6 7 8 Fig. 3. The expression profiles have been represented separately for each experimental group and patients arranged individually and C P, from treatment and control experiments, respectively. This relationship is defined as the significant intersection of genes between T P and C P, which is constrained by a threshold based on the typical statistical power of 80%. 2.4 Creation of Method Association Rule We create a set of method association rules that, given a set of differential profiles queried by the user, suggests the most appropriate method-associations capable to retrieve them. The method association rules are created based on the lattice structure from Section 2.2, containing all potential method-associations, and the set of all possible differential profiles P from Section 2.4 defined as })(,....,){( 1lCTCT PPPPP =where PPP jCT ∈)( represents each of the differential profiles present in P. 2.4.1 Method-Association Performance Evaluation We evaluate the performance of the method-associations MM i∈ for the query profiles ),,..,(1s SxxX = over two objectives: specificity and sensitivity )/( FNTPTNySpecificit += )/( FNTPTPySensitivit += , (2) where TP stands for True Positives (i.e., genes exhibiting profile S uXx ∈, which have been successfully retrieved by the applied method-association i M), TN stands for True Negatives (i.e., genes exhibiting profile S uXx ∉ and not retrieved by i M), FP stands for False Positives (i.e., genes exhibiting profile S uXx ∉and retrieved by i M) and FN stands for False Negatives (i.e., genes exhibiting profile S uXx ∈and not retrieved by i M). These four factors are calculated as: u iu TP ϕ ηϕ I =)(D )η(D)(D TN u iu ϕ ϕ − −− =I u iu u iu )η(D FN )(D η)(D FP ϕ ϕ ϕ ϕ − = − − =II , (3) where u ϕ represents the genes in the microarray set D that exhibit the queried profile S uXx ∈, and (D)Mii = η , the genes from D retrieved by the method-association i M. 2.4.2 Method-Association Selection We evaluate the method-associations in M based on their specificity and sensitivity. These two objectives are always conflicting, so we use a multiobjective optimization technique to maximize them, allowing us to detect all optimal methods-associations in M for the query profiles S X [11], [12]. We define objectives )( 2,1 OO corresponding to specificity and sensitivity respectively. 2.4.3 Creation of a Set of Method Association Rules We use the non-dominated method-associations described in Section 2.4.2 to create the method association rules R},...,{ 1k RR= where R R f∈is defined as: f R : IF 1 xIS f CT PP 1 )( AND , . . . , AND s x IS f sCT PP )( THEN f z IS i MWITH f C, (4) where )xx s ,...,( 1 are the profiles S X queried by the user; f CT PP 1 )(,…, PPP f sCT ∈)( ; Mz f∈ is the appropriate method-association to retrieve S X according to rule f R ; and f C denotes a measure of the specificity/sensitivity levels for f z, defined as: 21 2211 ))(*())(*( ww MOwMOw C ii f + + =, (5) where 1 w and 2 w are the weights associated to ),( 21 OO respectively. These values are provided by the user based on the relevance of each of these objectives for the particular study. If no values are given, the standard (0.5, 0.5) are used. 2.5 Prediction Using Method Association Rules The prediction phase works at two levels depending on the given input. If the input is a microarray data set D’, our methodology will provide the differential expression profiles P’ in the data set along with the optimal method-associations to retrieve such profiles. It might be the case that some of the differential profiles P’ uncovered from D’ were not included in the set of differential profiles P already learned by the methodology. Consequently, the information provided as input will be used to update P and R. If the input is a set of query profiles S X, the output will consist of the optimal method-association h M for S X at a certain f C value. To obtain these outputs, we apply matching and inference operations to the method association rule set [17]. Given an association rule set },...,{ 1k RRR =, for the differential profiles provided as the query set ),...,( 1s SxxX =, we define the matching degreeQof S uXx ∈with the if-part of the association rule f R as: f uCTu f uCTu PPxPPxQ )(1))(( ,−−= , (6) with being the Euclidean distance, and )( CT PP the centroids of the profiles. Therefore, given a set of query profiles S X, we define the strength of activation of the if-part of the rule f R as: )))(,(,...,))(,(min()( 1 1f sCTs f CT Sf PPxQPPxQXR =. (7) Let )),(( fSff CXRh denote the degree of association of the query profiles S X with the method-association i Maccording to rule f R and the specificity/sensitivity level f C. This degree is obtained by applying a product operator between )(Sf XR and f C. The optimal method-association for the queried profiles S X is defined as: i M/( i hkf iSi CXR ∈ =max)),(( f h)),(fSf CXR . (8) 3 Results We apply our procedure to a data set derived from longitudinal blood expression profiles of human volunteers treated with intravenous endotoxin compared to placebo. We expect to identify molecular pathways that provide insight into the host response over time to systemic inflammatory insults, as part of a Large-scale Collaborative Research Project sponsored by the National Institute of General Medical Sciences (www.gluegrant.org) [18]. The data were acquired from blood samples collected from eight normal human volunteers, four treated with intravenous endotoxin (i.e., patients 1 to 4) and four with placebo (i.e., patients 5 to 8) [18]. Complementary RNA was generated from circulating leukocytes at 0, 2, 4, 6, 9 and 24 hours after the i.v. infusion and hybridized with GeneChips® HG-U133A v2.0 from Affymetryx Inc., containing a set of 22283 genes. 3.1 Identification of Differentially Expressed Genes The statistical methods harbored have been applied using the standard p-value 05.0 =. The number of differentially expressed genes retrieved by each of the methods from the original set of genes is 1 M-10942 genes, 2 M-7841, 3 M-3904, 4 M-8023, 5 M-13151, 6 M-4588, 7 M-6070, 8 M-8557, 9 M-3995, 10 M-3367. These values show the number of significant genes retrieved by each of the statistical methods ranges in a wide rank. Moreover, the concordance rates also vary widely, inTable 1. Coincidence between methods in the retrieval of genes. The number in each cell represents a ratio of coincidence between genes retrieved by the statistical method in that column and the genes retrieved by the statistical method in that row relative to the total number of genes retrieved by the method in the row ( RowColumnRow /)( I). % 1 M2 M3 M 4 M 5 M6 M7 M 8 M9 M10 M 1 M -- 92.20 52.29 75.05 96.48 69.23 85.55 70.06 61.33 50.52 2 M56.06 -- 34.07 57.84 85.27 59.54 71.11 62.64 50.57 42.98 3 M 82.19 88.07 -- 96.24 94.77 57.35 78.75 72.87 56.86 46.73 4 M67.22 85.19 54.84 -- 95.16 55.49 73.65 70.20 51.49 42.83 5 M55.20 77.80 33.45 58.94 -- 50.28 66.72 66.38 46.42 38.93 6 M59.04 83.51 31.11 52.84 77.30 -- 89.63 56.56 60.64 49.38 7 M58.36 79.79 34.18 56.10 82.05 71.70 -- 62.34 57.23 49.07 8 M 57.36 84.34 37.96 64.17 95.96 54.30 74.80 -- 49.62 40.51 9 M62.10 84.21 36.63 58.21 84.74 72.00 84.95 61.36 -- 72.31 10 M 59.56 83.34 35.05 56.37 82.72 68.26 84.80 58.34 84.19 -- σ dicating that none of the methods subsumes the others (Table 1)(e.g., from the genes retrieved by 3 M, only 31.11% are also retrieved by 5 M, and 52.29% by 1 M). 3.2 Association of Statistical Methods The lattice arranged in this particular work contains all potential combinations of union and intersection of the ten statistical methods applied. Thus, M’ is defined as }...,...,,..., ,...,,,.,,..,{ 10932132132 31211021 MMMMMMMMMM MMMMMMMM ⊕⊕⊕⊕⊕⊕⊕⊕ ⊕⊕= We found that there is a relationship between the statistical methods and the differential profiles they are able to identify (see Section 2.2), having differential profiles identified by some methods and not by others. For example, the differential profile in (Fig. 4(a)) harbors 29 genes in our dataset D and is only retrieved by those statistical methods that take into account the time factor (e.g., 2 M, which retrieves more than 90% of these genes). This happens because the statistical methods that consider the treatment vs. control factor make an average of the expression values from patients 1 and 2 with those of patients 3 and 4 by considering them as replicas. Consequently, the differential behavior between them is lost. T rea tmen t - 8 Treatment - 21 Control - 9 Contro l - 7 TREATMENT CONTROL a) b ) Fig. 4. Examples of differential profiles only identified by some of the statistical methods 3.3 Identification of Differential Profiles The expression profiles have been represented separately for each experimental group (Section 2.3), and patients arranged individually. In our current problem, with eight patients, four treated with intravenous endotoxin (i.e., patients 1 to 4) and four with placebo (i.e., patients 5 to 8), and data retrieved over time at hours 0, 2, 4, 6, 9 and 24, each profile is represented by 24 consecutive time points (see Fig. 5). The differential profiles extracted from the treatment group show different levels of expression change. For example, there are sets of genes sharing very high variations in the levels of expression (e.g., profiles 15, 19, 21, and 22 in Fig. 5). In addition, some other profiles show differential characteristics for the patients (e.g., profiles 8 and 16 in Fig. 5). In the control group, the profiles are more homogeneous than in the treatment group. Typically, testing the coincidence among different data sources and clustering methods serves as a tool to investigate the validity of the identified groupings [19]. We follow this guideline to increase the confidence in the obtained differential profiles. Therefore, we calculate the coincidence between our retrieved differential profi- Treatment - 1 Treatment - 4 Treatment - 7 T re atme nt - 10 Treatment - 2 Treatment - 5 Treatment - 8 Treatment - 11 T re atme nt - 3 T re atme nt - 6 T re atme nt - 9 Treatment - 12 0 5000 10000 15000 20000 25000 0 5000 10000 15000 20000 25000 0 5000 10000 15000 20000 25000 0 5000 10000 15000 20000 25000 Tre atment - 13 Tre atment - 16 Tre atment - 19 Tre atment - 22 Tre atment - 14 Tre atment - 17 Tre atment - 20 Tre atment - 23 Tre atment - 15 Tre atment - 18 Tre atment - 21 Tre atment - 24 Co ntrol - 1 Co ntrol - 4 Co ntrol - 7 Control - 10 Co ntrol - 2 Co ntrol - 5 Co ntrol - 8 Control - 11 Control - 3 Control - 6 Control - 9 Control - 12 0 5000 10000 15000 20000 25000 0 5000 10000 15000 20000 25000 0 5000 10000 15000 20000 25000 0 5000 10000 15000 20000 25000 Treatment - 1 TREATMENT CONTROL Fig. 5. Representation of the differential profiles obtained separately for the treatment and control groups using the statistical methods applied in the current work les and external information provided by the Gene Ontology database [20]. To address this problem we developed an evolutionary multiobjective conceptual clustering methodology (R.R.Z., C.R.E., O.C., J.P.C., and I.Z., manuscript in preparation) that extracts clusters composed of features such as biological processes, molecular functions and cellular components defined at different specificity levels, and compare these clusters with our differential profiles by using a coincidence index test based on the hypergeometric distribution [9], [10], [19]. 3.4 Creation of Method Association Rules We have arbitrarily selected six profiles (i.e., 1 )( CT PP ,…, 6 )( CT PP ) identifying a total of 1395 genes in our dataset D and plotted as treatment clusters 2, 3, 4, 5, 10 and 12 in Fig. 5. These profiles represent genes exhibiting non-uniform behavior for distinct patients in the treatment group, and genes with changes in a level of expression smaller than 5000. We applied our methodology to find the optimal method-associations i M to retrieve them. 3.4.1 Method Association Performance Evaluation The results of the evaluation of the method-associations contained in the lattice M’ for the differential profiles are shown in Table 2, where the information relative to the sensitivity and specificity levels for the application of the most representative methodassociations over D is also specified. On the one hand, we observe that the union set of the genes obtained by seven of the statistical methods evaluated (i.e., methods 10876532 ,,,,,, MMMMMMM ) contains the 1395 genes desired (i.e., sensitivity value of 1) but with a low level of specificity (i.e., value of 0.369). On the other hand, the intersection set of genes obtained by the same seven statistical methods has a very low level of sensitivity (i.e., only 95 out of the 1395 genes were retrieved), whereas the value for specificity is very high. In between these two extremes we see some other method-associations which evaluation reveal trade-off solutions between the specificity and sensitivity objectives (Table 2). 3.4.2 Method Association Selection Once the method-associations M have been evaluated, we search for the nondominance relations in their applications to the microarray dataset D. The decision is based on the levels of specificity and sensitivity in Table 2. The Pareto optimal front conformed by this set of non-dominated method-associations is represented in Fig. 6. Table 2. Specificity and sensitivity values for the method-associations. The non-dominated solutions are pointed out with a star. Methods Specificity Sensitivity 2 M0.611 0.707 3 M0.826 0.205 5 M0.448 0.785 * 6 M0.813 0.447 * 7 M0.747 0.587 8 M0.625 0.537 * 10 M0.859 0.322 2 M∩3 M0.803 0.432 *2 M ∪ 3 M0.618 0.866 *Union of ( ),,,,,, 10876532 MMMMMMM 0.3690 1 *Intersection of ( ),,,,,, 10876532 MMMMMMM 0.983 0.066 1 0 0,2 0,4 0,6 0,8 1 00,2 0,4 0,6 0,8 1 Sensitivity Specificity (0.066, 0.983) (0.322, 0.859) (0.447, 0.813) (0.587, 0.747) (0.866, 0.618) (1, 0369) Fig. 6. Results of the evaluation of the method-associations contained in the lattice M’ for the six selected differential profiles 3.4.3 Creation of Method Association Rules The set of method association rules is created based on the evaluated profiles (i.e., 1 )( CT PP ,…, 6 )( CT PP ), and the method-associations i Mpresent in the Pareto optimal front of non-dominated solutions. The weights ),( 21 ww associated to the objectives ),( 21 OO are set to (0.5, 0.5) to calculate the specificity/sensitivity measure f C. We illustrate two association rules extracted from the evaluation of M’ over the former profiles, which have the following form: 1 R : IF 1 x IS 1 1 )( CT PP AND ,…, AND 6 x IS 1 6 )( CT PP THEN 1 Z IS 6 MWITH 1 C where f Cis calculated based on the specificity/sensitivity levels obtained on the application of such method over 1 )( CT PP ,…, 6 )( CT PP profiles (Table 2):