Prediction of Mitochondrial Matrix Protein Structures Based on Feature Selection and Fragment Assembly

Prediction of Mitochondrial Matrix Protein Structures Based on Feature Selection and Fragment Assembly Gualberto Asencio-Cort´es1,Jes´us S. Aguilar-Ruiz1, Alfonso E. M´arquez-Chamorro1,RobertoRuiz 1, and Cosme E. Santiesteban-Toca2 1School of Engineering, Pablo de Olavide University, Seville, Spain {guaasecor,aguilar,amarcha,robertoruiz}@upo.es 2Centro de Bioplantas, University of Ciego de ´ Avila, Cuba [email protected] Abstract. Protein structure prediction consists in determining the thre- e-dimensional conformation of a protein based only on its amino acid sequence. This is currently a difficult and significant challenge in structural bioinformatics because these structures are necessary for drug designing. This work proposes a method that reconstructs protein structures from protein fragments assembled according to their physico-chemical similarities, using information extracted from known protein structures. Our prediction system produces distance maps to represent protein structures, which provides more information than contact maps, which are predicted by many proposals in the literature. Most commonly used amino acid physico-chemical properties are hydrophobicity, polarity and charge. In our method, we performed a feature selection on the 544 properties of the AAindex repository, resulting in 16 properties which were used to predictions. We tested our proposal on 74 mitochondrial matrix proteins with a maximum sequence identity of 30% obtained from the Protein Data Bank. We achieved a recall of 0.80 and a precision of 0.79 with an 8-angstrom cut-off and a minimum sequence separation of 7 amino acids. Finally, we compared our system with other relevant proposal on the same benchmark and we achieved a recall improvement of 50.82%. Therefore, for the studied proteins, our method provides a notable improvement in terms of recall. Keywords: Protein structure prediction, physico-chemical amino acid properties, fragment assembly, protein distance map, feature selection. 1 Introduction Knowing the protein native 3D structures is currently a difficult and significant challenge because these structures determine protein function and they are necessary to design new drugs. Experimental methods to determine protein structures, generally X-ray crystallography and nuclear magnetic resonance, are very expensive and they have limitations with the structures of some proteins. Moreover, the great number of protein sequences whose three-dimensional structures must be determined, make computational methods of protein structure prediction a very useful tool. Protein structure prediction (PSP) consists in determining a three-dimensio- nal model based only on the amino acid sequence of a protein and it is currently an issue with great significance in structural bioinformatics [1]. There are currently two main approaches for the PSP problem. The first is the ab initio methods, which find the structure that corresponds to a global minimum of a function, generally a energy function, based in sequence properties. These methods do not use any protein as a template, their computational cost is generally very high and their reliability decreases when the sequence length increases [2]. The second main approach is homology methods, also known as comparative modeling, which try to solve the structure based on protein templates (templatebased modeling). This approach is based on the structural conservation of proteins in a protein family, since the 3D structures are more conserved in evolution than sequences. These methods are considered the most currently reliable approach for PSP problem [2]. Template-based modeling methods achieve good results when solved structures are available for proteins with sequences similar to the sequence of the target protein. However, when no homologous proteins with solved structures exist, free-modeling is used. Within free-modeling methods we find the fragment assembly methods that reconstruct the structure of a protein from structural fragments of other proteins. Three of most relevant fragment assembly-based methods are Fragment-HMM [3], FragFold [4] and ROSETTA [5]. ROSETTA uses a two-stage approach, which begins with a low-resolution model and continues with a representation of all the atoms of the protein, with the goal of minimizing the corresponding energy function. Since all information used in structure prediction must be inferred from amino acid sequence, there is many useful information derived from sequence used in the literature. Among this information, there are recent methods that use a great set of physico-chemical properties of amino acids [6]. However, the most commonly used properties are hydrophobicity, polarity and charge, which are used, for example, in the models HP and HPNX [7]. There is a database of amino acid properties named AAindex [8] which contains currently 544 properties, from which we selected a subset of 16 in this work by a feature selection process. The motivation for applying feature selection (FS) techniques has shifted from being optional to becoming a real prerequisite for model building. Specifically, in the PSP problem, the feature selection was also applied and improves the accuracy of predictions [9]. Theoretically, having more features should give us more discriminating power. However, this can cause several problems: increased computational complexity and cost; too many redundant or irrelevant features; and estimation degradation in the classification error. Based on the generation procedure, FS can be divided into individual feature ranking (FR) and feature subset selection (FSS) [10,11]. FR measures featureclass relevance, then ranks features by their scores and selects the top-ranked ones. These methods are widely used due to their simplicity, scalability, and good empirical success [12]. However, FR is criticized because it can only capture the relevance of features to the target concept, while redundancy and basic interactions between features are not revealed. Furthermore, it is difficult to determine the number of features retained, because a threshold is required. In contrast, FSS attempts to find a set of features that performs well. It integrates the metrics for measuring feature-class relevance and feature-feature interactions. In this work, a hybrid algorithm was used, BARS [13], in order to handle large datasets to take advantage of the above two approaches (FR, FSS) [14]. This method decouple relevance analysis and redundancy analysis, and have proven to be more effective than ranking methods and more efficient than subset evaluation methods in many traditional high-dimensional datasets. There are many PSP algorithms currently in the literature that produce a contact map to represent the predicted structure [6,15]. In contrast, our method produces a distance map, which includes more information than a contact map because it incorporates the distances between all of the amino acids in the molecule, irrespective of whether they make contact. There are fewer proposals in the literature that predict distance maps [16], because it is more difficult to perform regression than classification (continuous distances instead binary contacts). Some authors discretize the distances to predict, providing an intermediate representation between contacts and continuous distances, such as the proposal of Walsh et al. 2009 [2] which uses 4-class distance maps. However, unlike 3D models, both distance and contact maps have the desirable property of being insensitive to rotation or translation of the protein molecule. Our method is a free-modeling approach based on fragment assembly that selects the best distances between pairs of amino acids using fragments of known structures of proteins. These fragments are chosen through a searching process for nearest neighbors by similarity in 16 physico-chemical properties of amino acids selected from the AAindex repository. We tested our methodology by performing predictions on mitochondrial matrix proteins from the Protein Data Bank (PDB) [17] with a maximum sequence identity of 30%. We have performed predictions with a minimum sequence separation of 7 amino acids, as has been used in the literature [18]. Finally, we compared our system with RBFNN method proposed by Zhang et al. 2005 [19] with the same proteins in the same experimental conditions. In section 2, we define the elements, procedures and evaluation measures used by our prediction method. In section 3, we detail the used protein datasets, the experimental settings and the achieved results. Finally, in section 4, we describe the main conclusions of the performed study and we outline approaches for future studies. 2 Methods 2.1 Representation of Protein Structures The representation of protein structure that we used is the distance map, which is a square matrix of order L,whereLis the number of amino acids in the protein sequence. The distance matrix is divided in two parts: observed part (upper triangular) and predicted part (lower triangular). The element (i, j), where i<j, of the distance matrix is the actual distance measured in angstroms (˚ A) between the amino acids ith and jth in the sequence. To measure the distances between amino acids, it is necessary to use a reference atom of each amino acid. The most commonly used reference atoms are the alpha carbon and the beta carbon of amino acids [18]. In our method, we used the beta carbon (with the exception of glycine, for which the alpha carbon was used). The distances predicted by the algorithm are stored in the lower triangular of the distance map. Thus, the element (i, j)withi>jof the distance matrix is the predicted distance measured in angstroms between the amino acids ith and jth of the protein sequence. 2.2 Construction of Protein Fragments Knowledge Base Our prediction system ASPpred (Amino acid Subsequences Properties-based Predictor), works in two phases. In the first phase, it constructs a gallery of protein fragments from all the subsequences of all the proteins in the training set. In the second phase, the target structures of the proteins in test set were predicted using the generated protein fragments model. The knowledge base consists of a set of vectors called prediction vectors. Each one of these vectors was obtained from one training protein subsequence and contains the physico-chemical properties of the amino acids ends of such fragment. The vector also contains the actual distance between them. In order to define our prediction vectors formally, it is necessary to define the following elements. In first place, an amino acid sequence of length Lis defined by s1...s L. A fragment or subsequence into a sequence is represented by s1...s b...s e...s L,wheresb...s eis the fragment, sbis the beginning amino acid of the fragment, seis its ending amino acid and 1 ≤b<e≤L. Moreover, physico-chemical properties are defined by P1...P m,wheremis the number of properties used by the algorithm. The value of the property Pi of an amino acid sjis defined by Pi(sj). The prediction vector of a fragment is defined by the tuple showed in Equation 1. {B1,E 1,...,B m,E m,D}(1) Where Dis the distance between amino acids sband se.Biand Eiare defined in Equations 2 and 3, respectively. Birepresents the physico-chemical distribution of the entire sequence with decreasing weighting starting at the first amino acid of the fragment. Eiis analogous to Bistarting at the last amino acid of the fragment. Bi=Pi(sb)+ L  j=1 j=b Pi(sj) L|b−j|,∀i∈{1..m}(2) Ei=Pi(se)+ L  j=1 j=e Pi(sj) L|e−j|,∀i∈{1..m}(3) Note that prediction vectors represent fragments of different lengths, but these lengths is not included in them. The physico-chemical properties included in the prediction vectors are explained in the next subsection. From the point of view of data mining, Biand Eiare the attributes of training instances and Dthe class to predict. 2.3 Physico-chemical Feature Selection To the aim of using the smallest and most effective set of physico-chemical properties, we performed a feature selection from the repository AAindex of physico-chemical properties of amino acids. This repository currently contains 544 amino acid properties. We used BARS to perform the feature selection over all the properties in AAindex. BARS is an agglomerative algorithm due to the way it constructs the final subset of selected features. The method begins by generating a ranking. Then, pairs of features are obtained with the ranking’s first features, in combination with each one of the remaining features on the list. The pairs of features are ranked according to the value of the evaluation, and the process is repeated, that is, the subsets made up by the first sets on the new list are compared with the rest of the sets. At the end, the algorithm returns the best positioned feature subset of all the subsets evaluated. BARS can use any measure to evaluate feature subsets. Taken into account this domain with a numeric class attribute where distance between amino acids is represented, we used linear regression as evaluator criteria when the search process is carried out to find a relevant and not redundant subset of features. The dataset that we used for the feature selection is published by Fariselli et al. 2001 [18], that contains 173 proteins with a sequence identity lower than 25%, without chain breaks and with alignments with more than 15 sequences in the corresponding families. This process results on 16 physico-chemical properties that are showed in Table 1 with the same name and description used in AAindex. 2.4 Structure Reconstruction The second phase of our system consists in obtaining the prediction vectors of the target proteins and in performing a full sequential search to compare each test prediction vector with the training prediction vectors achieved in the first phase. The objective was to find the training prediction vector most similar to Table 1. The 16 physico-chemical properties of amino acids considered from AAindex CHOC760104 Proportion of residues 100% buried LEVM760104 Side chain torsion angle phi(AAAR) MEIH800103 Average side chain orientation angle PALJ810107 Normalized frequency of alpha-helix in all-alpha class QIAN880112 Weights for alpha-helix at the window position of 5 WOLS870101 Principal property value z1 ONEK900101 Delta G values for the peptides extrapolated to 0 M urea BLAM930101 Alpha helix propensity of position 44 in T4 lysozyme PARS000101 p-Values of mesophilic proteins based on the distributions of B values NADH010102 Hydropathy scale based on self-information values in the two-state model (9% accessibility) SUYM030101 Linker propensity index WOLR790101 Hydrophobicity index JACR890101 Weights from the IFH scale MIYS990103 Optimized relative partition energies - method B MIYS990104 Optimized relative partition energies - method C MIYS990105 Optimized relative partition energies - method D each test prediction vector. For the search process, we consider only training vectors with the same amino acid ends (first and last of each subsequence) than the test vectors. Figure 1 shows this search scheme. test Bt 1Et 1... Bt mEt m? training . . . Br 1Er 1... Br mEr mDr . . . Fig. 1. Search for the most similar training prediction vector In the search scheme of the Figure 1, Bt 1...Bt mand Et 1...Et mare the elements of the test subsequence explained above and Br 1...Br mand Er 1...Er mare those of the training subsequence with more similarity to the test subsequence. The distance field Drof the most similar training vector is assigned to the distance field (symbolized with ?) of test vector. The training vector with the highest similarity to a test vector satisfies the condition showed in Equation 4. As can be seen in that condition, for the comparison of the test and training vectors, an Euclidean distance is used, which includes all the attributes in these vectors with same weigths. All these attributes are normalised previously. The normalization ensured that all of the attributes are on the same scale and contributed equally to the prediction. min r∈T rainingSet     m  i=1 (Bt i−Br i)2+ m  i=1 (Et i−Er i)2(4) Finally, once predicted distances are assigned in test vectors, these distances are stored in the lower triangular of the distance map of the test sequence. Specifically, the distance field of the prediction vector of the subsequence sb...s e is assigned to the position (e, b) of the distance map. Thus the structure of each target sequence, by its distance map, is reconstructed. 2.5 Evaluation of Predicted Models We used several measures to evaluate the quality of the predictions. The first measure is the precision, that is the percentage of predicted contacts that are present in the native structure. This measure is largely used in the literature of protein structure prediction, as in the works of Fariselli et al. [18,20]. The second one is the recall, that is the percentage of native contacts that are predicted to be contacts. Recall has also been widely used in other protein prediction methods [19]. Finally, we have obtained measures of accuracy, specificity and Matthews Correlation Coefficient, that may often provide a much more balanced evaluation of the prediction than percentages [21]. The following formulas (5,6,7,8,9) define these five measures. Precision =TP TP +FP (5) Recall =TP TP +FN (6) Accuracy =TP +TN TP +FP +FN +TN (7) Specificity =TN TN +FP (8) MCC =TP ×TN −FP ×FN (TP +FP)(TP +FN)(TN +FP)(TN +FN)(9) These measures are used to evaluate the quality of a classification: i.e., each predicted value is assigned a value of 0 or 1. Thus, there are four possible outcomes depending on the quality of the predictions: a) both the real and predicted values are 1 (true positive, TP), b) both the real and predicted values are 0 (true negative, TN), c) the real value is 1 and the predicted value is 0 (false negative, FN) and d) the real value is 0 and the predicted value is 1 (false positive, FP). Because in this case, the class to predict is a real value (a distance), to obtain these measures it is necessary to binarise the class using a distance threshold or cut-off. In this work, we used a cut-off value of 8 angstroms, which is commonly used in the literature [18,20,19]. In the evaluation of the measures, we omitted predictions of amino acid pairs with a minimum separation in the protein sequence of 7 amino acids, as in Fariselli et al. [18]. 3 Experimentation and Results 3.1 Prediction of Mitochondrial Matrix Proteins We performed an experimental validation of our predictor using all mitochondrial matrix proteins (GO ID: 5759) published in the PDB with a maximum identity of 30% (non-homologous proteins) at October 2011 (74 proteins with a maximum length of 1094 amino acids). In Table 2, we show the PDB codes of the proteins used in this study. We classified proteins in three groups of sequence length L (L≤300, 300 <L≤450 and L>450) in order to show the prediction behavior for each sequence length interval. Table 2. Mitochondrial matrix proteins used to train and test our predictor L≤300 2CW6 1CSH 2DFD 3CMQ 1PJ3 3C5E 1BWY 2GRA 1D2E 2E0A 3EXE 1WDK 3DLX 1EFV 2HDH 1F0Y 2IB8 3GH0 1WLE 3E04 1KKC 2O23 1GKZ 2IZZ 3KGW 1ZMD 3IHJ 1MJ3 2WYT 1HW4 2OAT 7AAT 2FGE 3IKL 1QQ2 3ED7 1I4W 2QB7 L>450 2J6L 3IKM 1R4W 3EMN 1OTH 2QFY 1A4E 2JDI 3MW9 1RHS 3QUW 1RX0 2R2N 1CJC 2UXW 3N9Y 1TG6 3ULL 1W6U 3AFO 1G5H 2WYA 3OEE 1XX4 5CYT 2A1H 3BLX 1HR6 2XIJ 3OU5 1ZD8 L300 −450 2BFD 3BPT 1OHV 2ZT5 3SPA A cross-validation was performed over each group of proteins and over the all 74 proteins. We used a leave-one-out scheme in order to avoid the effect of fold choice in a cross-validation with folds. Table 3 shows the five evaluation measures obtained in this experiment. We achieved a recall value of 0.80 and a precision of 0.79 for the complete group of proteins, as shown in Table 3. We obtain best predictions, in terms of recall and precision, with proteins of length between 300 and 450 amino acids. In this group of proteins, we achieved recall of 0.84 and precision of 0.83. In most cases the precision obtained in predicting proteins with long sequences (more than 300 amino acids) is lower than with proteins of short sequences. For example, in the work of Fariselli et al. 2001 [18], which also uses a cross validation, cut-off of 8 angstroms and minimum sequence separation of 7 amino acids, achieved a precision value of 0.11 for proteins of more than 300 amino acids. Table 3. Efficiency of our method predicting mitochondrial matrix proteins Protein set Recall Precision Accuracy Specificity MCC All proteins (74) 0.80 0.79 0.97 0.97 0.82 L≤300 (20) 0.77 0.76 0.98 0.98 0.75 300 <L≤450 (27) 0.84 0.83 0.99 0.99 0.83 L>450 (27) 0.77 0.76 0.95 0.95 0.82 (a) 1TG6 (277 amino acids) (b) 3BLX (349 amino acids) (c) Color scale Fig. 2. Predicted distance maps for the mitochondrial matrix proteins 1TG6 (a) and 3BLX (b) with their color scale (c) Figure 2 shows the predicted distance maps for protein 1TG6 (277 amino acids) and 3BLX (349 amino acids) from the study set. We show a color scale to represent the distances, ranging from the minimum (red) to the maximum (blue) distance. We can appreciate in these distance matrices that the lower triangular (predicted distances) is largely similar to the upper triangular (real distances). 3.2 Comparison with RBFNN on the Same Benchmark In order to assess the quality of the predictions obtained with our method and to validate our predictor, we compared our proposal with RBFNN method proposed by Zhang et al. 2005 [19]. We predicted the same test proteins with the same training sets in the same conditions. Zhang et al. 2005 used recall (namely accuracy (Ap) by the authors), predicted and desired numbers to evaluate the performance. Predicted numbers Npis the count of the predicted contacts by the algorithm and desired numbers Ndis the total number of contacts. The contact threshold was set at 8 ˚ A. In Table 4 we show the results of this experimentation. As we can see in this table, the average recall (Ap) of ASPpred is 50.82% higher than RBFNN.