Optimization of multi-classifiers for computational biology: application to gene finding and expression

See discussions, stats, and author profiles for this publication at: https://www.researchgate.net/publication/225450776 Optimization of Multi-classifiers for Computational Biology: Application to the Gene Finding Problem ArticleinTheoretical Chemistry Accounts · March 2009 DOI: 10.1007/s00214-009-0648-3·Source: IEEE Xplore CITATIONS 11 READS 99 4 authors: Some of the authors of this publication are also working on these related projects: Gene expression kinetics View project Improving Breath Test sensitivity through Neural Networks and other Machine Learning systems. View project Rocío Romero-Zaliz University of Granada 53 PUBLICATIONS753 CITATIONS SEE PROFILE Cristina Rubio-Escudero Universidad de Sevilla 70 PUBLICATIONS564 CITATIONS SEE PROFILE Igor Zwir University of Granada 135 PUBLICATIONS1,844 CITATIONS SEE PROFILE Coral del Val University of Granada 89 PUBLICATIONS2,077 CITATIONS SEE PROFILE All content following this page was uploaded by Coral del Val on 01 June 2014. The user has requested enhancement of the downloaded file. REGULAR ARTICLE Optimization of multi-classifiers for computational biology: application to gene finding and expression Rocı ´o Romero-Zaliz •Cristina Rubio-Escudero • Igor Zwir •Coral del Val Received: 17 March 2009 / Accepted: 23 September 2009 / Published online: 15 October 2009 The Author(s) 2009. This article is published with open access at Springerlink.com Abstract Genomes of many organisms have been sequenced over the last few years. However, transforming such raw sequence data into knowledge remains a hard task. A great number of prediction programs have been developed to address part of this problem: the location of genes along a genome and their expression. We propose a multi-objective methodology to combine state-of-the-art algorithms into an aggregation scheme in order to obtain optimal methods’ aggregations. The results obtained show a major improvement in sensitivity when our methodology is compared to the performance of individual methods for gene finding and gene expression problems. The methodology proposed here is an automatic method generator, and a step forward to exploit all already existing methods, by providing alternative optimal methods’ aggregations to answer concrete queries for a certain biological problem with a maximized accuracy of the prediction. As more approaches are integrated for each of the presented problems, de novo accuracy can be expected to improve further. Keywords Multiobjective Gene finding  Gene expression 1 Introduction Genomes of many organisms have been sequenced over the last few years. However, transforming such raw sequence data into knowledge remains a hard task [1]. A great number of prediction programs have been developed to address one part of this problem: the location of genes along a genome [2–4]. Unfortunately, finding genes in a genomic sequence is far from being a trivial problem. Computational gene prediction methods yet have to achieve perfect accuracy, even in the relatively simple prokaryotic genomes [1]. Gene prediction is one of the most important problems in computational biology due to the inherent value of the set of protein-coding genes for other analyses. Another part of the problem is determining when, where and for how long these genes are turned on or off. Microarray technology allows the simultaneous evaluation of the expression of hundreds of genes in a single assay, converting this technology in a powerful tool for expression profiling, as well as, diagnosis and classification of cancers and other diseases. However, this technology presents a wealth of analysis problems [5] such as the inherent variability of cDNA microarrays at the individual slide and spot level, the large-scale nature of the data, and the fact that the full use of expression profiles for inferring gene function is still only partly explored. Many new methods have been developed to address the statistical challenge of Dedicated to Professor Sandor Suhai on the occasion of his 65th birthday and published as part of the Suhai Festschrift Issue. Rocı ´o Romero-Zaliz, Cristina Rubio-Escudero contributed equally. R. Romero-Zaliz I. Zwir C. del Val (&) Computer Science and Artificial Intelligence Department, Granada University, Granada, Spain e-mail: [email protected] C. Rubio-Escudero Departamento de Lenguajes y Sistemas Informa ´ticos, Sevilla University, Sevilla, Spain I. Zwir Howard Hughes Medical Institute, Department of Molecular Microbiology, Washington University School of Medicine, St. Louis, MO, USA 123 Theor Chem Acc (2010) 125:599–611 DOI 10.1007/s00214-009-0648-3 identifying ‘‘important’’ genes in the large sets of raw sequence data [6–11]. However, there is still 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 of the many available statistical methods is the most reliable to identify differences across profiles. Despite the advances in both referred problems, existing approaches to predict genes and to analyze microarray data have intrinsic advantages and limitations [1,12]. Furthermore, there is no program or methodology that can provide perfect predictions for any given input for either of these two problems. The problems of gene finding (identifying genes, exons and introns, beginning and end of the genes) and analysis of gene expression are formulated in this paper as classification problems. The gene finding problem can be interpreted as a simple decision between which section of a sequence is protein coding and which is not. Concerning the gene expression from microarray experiments, the classification problem can be seen as a decision between which genes are active or inactive in a given time point and/or under a given condition. For both problems many different programs are available, which give distinct solutions. There have been previous approaches to combine gene predictors [13–15] and microarray analyses [16,17], but maximizing accuracy by weighting both sensitivity and specificity functions into a single objective. However, our methodology uses a multi-objective approach to extract the best methods’ aggregations by maximizing the specificity and sensitivity of their predictions individually. This approach combines state-of-the-art algorithms into an aggregation scheme to provide better predictions by taking advantage of the different methodologies’ strengths and avoiding their weaknesses. We applied our methodology to the both referred problems. In the gene finding problem, we used the EGASP sets from the ENCODE Genome Annotation Assessment Project (EGASP) [18,19]. These datasets contain manually curated fragments of the human genome originating from the ENCODE project [20]. This data set was selected by the EGASP assessment because the genes encoded in these regions were not used to train any particular gene predictor. Therefore, it is not a biased dataset. In the case of analysis of the microarrays, we used a dataset derived from the analysis of longitudinal blood expression profiles of human volunteers treated with intravenous endotoxin, compared to those treated with a placebo in order to study the inflammation and human response to injury. This dataset was part of a Large-scale Collaborative Research Project sponsored by the National Institute of General Medical Sciences [21]. 2 Materials and methods The aggregation of the results from various methods is accomplished using the union -[- and intersection -\- operators [22]. All potential aggregations, termed methods’ aggregations from here on, form a space of potential hypotheses, which can be represented as a lattice structure (Fig. 1). We search for the best methods’ aggregations, moving from hypothesis to hypothesis towards the most general, the union of all methods, and the most specific, their intersection, which are located at the top and the bottom of the lattice, respectively [23] (Fig. 1). In Fig. 1 we can appreciate the lattice generated by the union and intersection of methods. In the gene finding problem we explore five methods, n=5, termed M 1 to M 5 , forming a total set of 31 potential aggregations. In the analysis of microarrays problem, ten methods are surveyed, n=10, termed M 1 to M 10 , forming a total set of 1,023 potential aggregations. The methods’ aggregations are evaluated based on a multi-objective approach [24] to extract the best methods’ aggregations by maximizing the specificity (Sp) and sensitivity (Sn) of their predictions. To estimate the sampling bias [25] of the methods’ aggregations we randomly partition the original sample into ten subsamples, and each subsample was retained as a validation data for each methods’ aggregations. 2.1 Gene finding problem: dataset and programs For the gene finding problem, we selected 27 ENCODE regions to test our proposal. These ENCODE regions have undergone an exhaustive annotation strategy prior to EGASP by the HAVANA team [26]. They consist of 2,471 total transcripts representing 434 unique protein-coding gene loci. The programs used in this study are those used in the EGASP competition, which are ab initio gene predictors using a single genome sequence. These programs were designed to predict gene structure, or at least a set of spliceable exons in vertebrate or pre-human genome sequences: GeneID [27], Genscan [28], Genemark [29], Augustus [30] and GeneZilla [31]. GeneID combines different algorithms using Position Weight Arrays to detect features such as splice sites, start and stop codons and Markov Models to score exons and Dynamic Programming (DP) to assemble the gene structure [27]. Genescan uses a general probabilistic model for the gene structure of human genomic sequences. It has the capacity to predict multiple genes in a sequence, to deal with partial as well as complete genes, and to predict consistent sets of genes occurring on either or both DNA strands [32]. Genemark for eukaryotes gathers the original Genemark models into the naturally derived hidden Markov model framework with 600 Theor Chem Acc (2010) 125:599–611 123 gene boundaries modeled as transitions between hidden states [29]. Augustus is a gene predictor for eukaryotic genomic sequences that is based on a generalized hidden Markov model, a probabilistic model of a sequence and its gene structure [30]. GeneZilla is based on the Generalized Hidden Markov Model (GHMM) framework, similar to Genscan. Graph-theoretic representations of the high scoring open reading frames are provided, allowing for exploration of sub-optimal gene models. It makes use of Interpolated Markov Models (IMMs), Maximal Dependence Decomposition (MDD), and includes states for signal peptides, branch points, TATA boxes and CAP sites [31]. For each method, the closest organism available for each gene in the dataset was selected. Predictions on both strands were extracted. The aggregations of the results of the different gene prediction approaches are performed at a nucleotide level. The aggregation of the results of different methods joins two or more overlapping or adjacent exons into a larger new exon (Fig. 2). Nucleotide level accuracy is calculated as a comparison of the annotated nucleotides with the predicted nucleotides. Individual nucleotides appearing in more than one transcript in either the annotation or the predictions are considered only once for the nucleotide level statistics. Nucleotide predictions must be on the same strand as the annotations to be counted as correct. At the nucleotide level, Sn is the proportion of annotated nucleotides (as being coding or part of an mRNA molecule) that is correctly predicted, and Sp the proportion of predicted nucleotides (as being coding or part of an mRNA molecule) that is so annotated. As a summary measure, we have computed the correlation coefficient between the annotated and the predicted nucleotides [19,33]. Sp ¼TP TP þFP Sn ¼TP TP þFN ð1Þ CC ¼ðTP TNÞðFN FPÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðTP þFNÞðTN þFPÞðTP þFPÞðTN þFNÞ pð2Þ 2.2 Gene expression profile finding: datasets and analysis methods The dataset used was derived from longitudinal blood expression profiles of human volunteers treated with intravenous endotoxin compared to those treated with a placebo. The data are related to the host response over time to systemic inflammatory insults, as part of a large-scale collaborative research project sponsored by the National Institute Fig. 1 Lattice of potential hypothesis, methods’ aggregations of M 1 …M n using the -[- and -\- and operators. The solid arrows show the direction of the search in the space of hypotheses Theor Chem Acc (2010) 125:599–611 601 123 of General Medical Sciences (http://www.gluegrant.org). The data was derived from blood samples collected from eight normal human volunteers, four treated with intravenous endotoxin (i.e., patients 1–4) and four with placebo (i.e., patients 5–8) [21]. Complementary RNA was generated from circulating leukocytes at 0, 2, 4, 6, 9 and 24 h after the intravenous infusion and hybridized with GeneChips  HG-U133A v2.0 from Affymetryx Inc., containing a set of 22,283 probe sets. A total set of 29 gene expression profiles (sets of genes which exhibit a common behavior throughout the conditions of the problem under study, time, treatment and patient in our particular case) are contained in the dataset and forms the focus of our study [34]. The methods analyzed in this study were applied to identify meaningful gene expression profiles from microarray data. The list of programs used comprises the methods most frequently applied to the analysis of microarray data: Student’s ttest [7], Permutation Test [10], Analysis of Variance (ANOVA) [7] and Repeated Measures Analysis of Variance (RMANOVA) [35]. These methods have been applied to the inflammation and host response to injury problem to account for different experimental conditions, such as treatment versus control and different time points. Therefore, Student’s ttest and Permutation Test have been applied in two different ways: considering treatment versus control and considering time. The ANOVA and RMANOVA tests can account for more than one experimental condition simultaneously; therefore they have been applied in three different ways: considering treatment versus control, considering time, and considering treatment versus control and time simultaneously. The aggregation of the results of different methods in the Gene Expression Profile Finding Problem is performed by combining the results obtained by each of the individual methods (group of probe sets). The union of two methods, M a and M b (M a [M b ), is defined as the group resulting which contains all genes retrieved either by method M a or by M b. The intersection of two methods, M a and M b (M a \ M b ), is defined as the group resulting which contains all genes retrieved by both methods M a and M b. We evaluate the performance of each method aggregation to retrieve each of the 29 gene expression profiles present in the inflammation and host response to injury dataset. In our particular problem, when studying the behavior of method M i to retrieve a gene expression profile P j , we define true positives (TPs) as probe sets retrieved by method M i which exhibit the gene expression profile P j , true negatives (TNs) probe sets not retrieved by method M i which do not exhibit the gene expression profile P j , false positives (FPs) as probe sets retrieved by method M i which do not exhibit the gene expression profile P j , and false negatives (FNs) as probe sets not retrieved by method M i which exhibit the gene expression profile P j . TP, TN, FP and FN information is typically summarized in terms of Sn, the proportion of probe sets belonging to P j in the dataset and correctly retrieved by the method M i under evaluation, and Sp, the proportion of probe sets correctly retrieved by the method M i from all the probe sets Fig. 2 Aggregation of exons using the union operator aequal exons, bmissed exons, cincluded exons, doverlapped exons 602 Theor Chem Acc (2010) 125:599–611 123 retrieved by method M i (see Eq. 1). These measures are formally described as for the Gene Finding Problem. 3 Results The results obtained applying our methodology to the two proposed biological problems outperform in terms of specificity and sensitivity the results obtained by classical methods, though both gene prediction and the identification of gene expression profiles are problems of different nature. 3.1 Gene finding The most updated version of the ENCODE annotations and of each gene finding prediction algorithm was used. The specificity, sensitivity and correlation coefficient (CC) averages for each individual method are shown in Table 1 with values represented in the [0–1] interval (see Eqs. 1,2). Genscan showed the highest CC while GeneMark obtained the lowest CC. GeneID obtained the highest specificity and the lowest sensitivity, while GeneZilla showed the highest sensitivity with lower specificity. The analysis of the individual results shows that some algorithms are able to predict certain genes very accurately with CC values close to 1, but the same algorithm completely fails to predict other genes (CC below 0.7 or even 0.5). These results show that a high average CC does not imply a good performance, and vice versa, since the average might hide some low CCs for specific genes. We evaluated all possible 31 methods’ aggregations resulting from applying the union and intersection operators to the selected five gene finding programs. The results (Fig. 3) show the general increase of sensitivity and the decrease of specificity with the increasing number of methods per aggregation when applying the union operator. However, we find the opposite behavior when using the intersection operator. The comparison between the prediction accuracy of the individual methods and the aggregation strategy can be seen in Fig. 3. This figure shows that methods’ aggregations increase either prediction’s sensitivity or specificity, and a few aggregations increase both objectives when compared with the individual methods. The best results are obtained in aggregations containing 2 or 3 methods via either union or intersection. Figure 4shows the Pareto optimal front [36] in red, which is the best methods’ aggregation in terms of sensitivity and specificity simultaneously. The Pareto optimal front consists of those methods’ aggregations for which improvement in specificity can only occur with the worsening of sensitivity and vice versa, that is, the best method aggregations in terms of sensitivity and specificity, simultaneously. Individual methods are represented in blue. The top ten methods’ aggregations are shown in Table 2. The results show that aggregations improving the individual methods’ performances in sensitivity and specificity generally include Augustus. There are several best methods’ aggregations, including union and intersection, but the one requiring the lower number of methods and predicting the highest number of genes is the union Augustus [GeneID. The best combination is Genscan \ GeneZilla when considering the intersection operator. The former correctly predicts 80.00% of the dataset, while the latter achieves 73.33%. GeneZilla is present in many of the best ten methods’ aggregations, providing supplementary predictions to the other methods, even though it is the worst individual method (Table 1). Table 1 Individual gene finding method’s performance Method Specificity Sensitivity Correlation coefficient Genscan 0.772 0.759 0.745 Genzilla 0.750 0.782 0.744 Augustus 0.794 0.694 0.724 GeneID 0.829 0.658 0.712 GeneMark 0.764 0.664 0.683 The methods are ordered according to their correlation coefficient; the best result for each column is highlighted in bold and the worst in italic Fig. 3 Sensitivity versus specificity for all methods’ aggregations and individual methods. The number of methods used in each aggregation is color-coded as indicated in the legend. The five lightyellow circles represent sensitivity and specificity values for the individual methods Theor Chem Acc (2010) 125:599–611 603 123 Cross-validation techniques are often used to estimate how accurately a predictive model will perform in practice, and consequently, to specifically avoid overfitting [23]. Unfortunately, because the original training datasets and some programmed code of the individual methods are not accessible, they cannot be re-trained based on different data partitions and tested with the remaining ones. Therefore, to indirectly estimate the sampling bias of the methods’ aggregations [25], we randomly partition the original sample into ten subsamples, and each subsample was used as validation dataset for each methods’ aggregation (Table 2). The results obtained show a small sampling variability (Fig. 5), suggesting good performance and robust results for the methods’ aggregations selected as the best ones. The levels of specificity and sensitivity obtained by each methods’ aggregation over the complete dataset (ENCODE) are shown in Fig. 6, ranging from 0 (green) to 1 (red). Each row represents a methods’ aggregation and each column a gene from the dataset. The rows and columns for each graph were clustered independently, and therefore we can see groups (clusters) of method’s aggregations showing similar behavior. For instance, in Fig. 6a, there are several methods’ aggregations using the union operator that cannot predict some ENCODE genes; those are represented as green columns (e.g., AC068580, AC079630, AC021607). Regions, which are easy to detect by many of the methods’ aggregations are represented as mostly red columns (e.g., AC093511, AC131574, AC113331). The sensitivity of each method aggregation using the union operator to predict each gene is shown in Fig. 6b. There are a few genes (e.g., AC068580, AC079630, AC021607) that obtain very low sensitivity for all methods’ aggregations, as it can be seen in their green columns. The aggregation of methods increases the sensitivity of the prediction as it is shown by red cells in Fig. 6b. Figure 6c shows gene prediction specificity for each methods’ aggregation using the intersection operator. The methods’ aggregations using the intersection operator increase the specificity of the prediction when compared to single methods and the methods’ aggregations using the union operator, as it is illustrated in Fig. 6c. However, there are several genes that are not predicted by any method (e.g., AC068580, AC079630, AC021607). Sensitivity, on the other hand, has a different behavior for each methods’ aggregation (Fig. 6d). Some genes are more Fig. 4 Sensitivity versus specificity: Optimal Pareto front. There are three color-codes: blue dots correspond to single methods, red dots to those methods’ aggregations that are optimal in both objectives sensitivity and specificity (Pareto optimal front); and black dots are all other methods’ aggregations Table 2 Ten best methods’ aggregations Methods’ aggregations Specificity Sensitivity Correlation coefficient % Genes correctly predicted Augustus [GeneID 0.805 0.807 0.791 80.00 Augustus [Genscan 0.771 0.839 0.786 74.17 Genscan \GeneZilla 0.850 0.683 0.745 73.33 Augustus [Genscan [GeneID 0.757 0.851 0.782 71.67 Augustus [Genscan [GeneID [GeneMark 0.725 0.864 0.767 70.83 Augustus [GeneZilla [GeneMark 0.706 0.874 0.760 70.00 Augustus [Genzilla [Genscan [GeneID 0.700 0.882 0.756 69.17 Genzilla \GeneID 0.866 0.612 0.706 68.33 Genzilla \Genscan \GeneID 0.711 0.870 0.758 67.50 Augustus \Genzilla 0.851 0.640 0.714 66.67 Gene finding methods are ordered by descending number of genes correctly retrieved. A gene is considered correctly retrieved when its correlation coefficient is over 0.7. The best result for each column is highlighted in bold 604 Theor Chem Acc (2010) 125:599–611 123 difficult to predict than others as represented by mostly red (e.g., AC072051, AL023881) or green columns (e.g., AC079630, AC068580). Finally, there are several genes which are not recognized by any method or methods’ aggregations (e.g., AC068580, AC079630, AC021607). 3.2 Gene expression profile finding The profiles conforming the dataset are the ones obtained from the inflammation and host response to injury problem. The methods analyzed in this study are: Student’s ttest Fig. 5 Sub-sample boxplots. aEach boxplot represents the specificity obtained over all ten sub-sample sets from the original dataset with a method aggregation. bEach boxplot represents the sensitivity obtained over all sub-sample sets of the original dataset with a method aggregation. The methods’ aggregations applied are the ones reported in Table 2as the ten best methods’ aggregation Fig. 6 Specificity and sensitivity of methods’ aggregations obtained for the ENCODE dataset. Each column represents a gene of the ENCODE region (120 columns) and each row a method aggregation including the individual ones (31 rows). a Specificity for the union operator. bSensitivity for the union operator. cSpecificity for the intersection operator. d Sensitivity for the intersection operator. The color is coded from 0 (green)to1(red). Green cells represent low values while the red cells correspond to high values. Black points correspond to values around 0.5. Labels for the xand yaxis are shown in additional Table 1 Theor Chem Acc (2010) 125:599–611 605 123 considering treatment versus control, Student’s ttest considering time, Permutation Test considering treatment versus control, Permutation Test considering time, ANOVA considering treatment versus control, ANOVA considering time, ANOVA considering treatment and time, RMANOVA considering treatment, RMANOVA over time, and RMANOVA over treatment and time. These methods have been relabelled for simplification purposes (Table 3). We show in Table 4the average results to retrieve each of the 29 gene expression profiles obtained by the individual methods in terms of specificity, sensitivity and correlation coefficient (CC). Values are represented in the [0–1] interval. Out of all gene expression methods analyzed, ANOVA considering time (represented by M 5 ) achieved the best sensitivity level and the Permutation Test considering time (M 3 ), the best specificity level. However, their average correlation coefficient (CC) for all profiles was not the highest. We can see that the levels of CC are generally low. This is due to the type of problem we are dealing with, finding a particular profile in a very large set of data, and obtaining as a result a large rate of false positives (FP), which decrease dramatically the correlation coefficient associated to each method. However, as occurred with the gene prediction problem, the results show that a low average CC does not imply a bad performance. Some methods recover particular profiles with better levels of CCs, and also with a good specificity/ sensitivity aggregation levels. All 1,023 potential methods’ aggregations of the ten applied gene expression programs were evaluated. Figure 7 represents the general increase of sensitivity with the number of methods combined using the union operator, while specificity increases when combining with the intersection operator. Individual methods (in red), are in a middle position generally improved in terms of specificity/ sensitivity by methods’ aggregations. These methods’ aggregations provide a wide spectrum of specificity/sensitivity levels depending on the operator applied. Union favors sensitivity, decreasing the number of false positives (FP), whereas intersection favors specificity, decreasing the number of false negatives (FN). The correlation coefficient is improved by application of any of the two operators. We show in Table 5the ten best methods’ aggregations in terms of CC. We see how the CC levels widely overcome the levels obtained by the single methods. Methods aggregated by means of the union operator also highly improve the sensitivity levels, while methods obtained applying the intersection operator widely overcome the single ones in terms of specificity. The CC levels obtained Table 3 Relabelling of methods analyzed in this study Legend Methods M 1 Student’s ttest considering treatment versus control M 2 Student’s ttest considering time M 3 Permutation test considering treatment versus control M 4 Permutation test considering time M 5 ANOVA considering treatment versus control M 6 ANOVA considering time M 7 ANOVA considering treatment and time M 8 RMANOVA considering treatment M 9 RMANOVA over time M 10 RMANOVA over treatment and time Table 4 Results obtained by the individual methods, M 1 to M 10 Methods Specificity Sensitivity Correlation Coefficient M 1 0.553 0.560 0.019 M 2 0.266 0.776 0.007 M 3 0.716 0.342 0.011 M 4 0.502 0.551 0.009 M 5 0.195 0.818 0.003 M 6 0.477 0.530 0.002 M 7 0.346 0.656 0.001 M 8 0.455 0.546 0.001 M 9 0.559 0.451 0.002 M 10 0.621 0.392 0.002 The best result for each column is highlighted in bold, while the worst result is highlighted in italic Fig. 7 Sensitivity versus specificity for all methods’ aggregations and individual methods. The number of methods used in each aggregation is color-coded as indicated in the legend. The lightyellow circles represent sensitivity and specificity values of the individual methods 606 Theor Chem Acc (2010) 125:599–611 123