Improved residue contact prediction using support vector machines and a large feature set

Background Predicting protein residue-residue contacts is an important 2D prediction task. It is useful for ab initio structure prediction and understanding protein folding. In spite of steady progress over the past decade, contact prediction remains still largely unsolved. Results Here we develop a new contact map predictor (SVMcon) that uses support vector machines to predict medium- and long-range contacts. SVMcon integrates profiles, secondary structure, relative solvent accessibility, contact potentials, and other useful features. On the same test data set, SVMcon's accuracy is 4% higher than the latest version of the CMAPpro contact map predictor. SVMcon recently participated in the seventh edition of the Critical Assessment of Techniques for Protein Structure Prediction (CASP7) experiment and was evaluated along with seven other contact map predictors. SVMcon was ranked as one of the top predictors, yielding the second best coverage and accuracy for contacts with sequence separation >= 12 on 13 de novo domains. Conclusion We describe SVMcon, a new contact map predictor that uses SVMs and a large set of informative features. SVMcon yields good performance on medium- to long-range contact predictions and can be modularly incorporated into a structure prediction pipeline.


Background
Predicting protein inter-residue contacts is an important 2D structure prediction problem [1]. Contact prediction can help improve analogous fold recognition [2,3] and ab initio 3D structure prediction [4]. Several algorithms for reconstructing 3D structure from contacts have been developed in both the structure prediction and determination (NMR) literature [5][6][7][8]. Contact map prediction is also useful for inferring protein folding rates and pathways [9,10].
Due to its importance, contact prediction has received considerable attention over the last decade. For instance, contact prediction methods have been evaluated in the fifth, sixth, and seventh editions of the Critical Assessment of Techniques for Protein Structure Prediction (CASP) experiment [11][12][13][14][15]. A number of different methods for predicting contacts have been developed. These methods can be classified roughly into two non-exclusive categories: (1) statistical correlated mutations approaches [16][17][18][19][20][21][22]; and (2) machine learning approaches [23][24][25][26][27][28][29][30][31][32][33][34]. The former uses correlated mutations of residues to predict contacts. The latter uses machine learning methods such as neural networks, self organizing map, hidden Markov models, and support vector machines to predict 2D contacts from the primary sequence, as well as other 1D features such as relative solvent accessibility and secondary structure.
In spite of steady progress, contact map prediction remains however a largely unsolved challenge. Here we describe a method that uses support vector machines together with a large set of informative features to improve contact map prediction. On the same data set, SVMcon outperforms the latest version of the CMAPpro contact map predictor [28,35] and is ranked as one of the top predictors in the blind and independent CASP7 experiment.

Results and Discussion
We first compare SVMcon with the latest version of CMAPpro on the same benchmark dataset. Then we describe the performance of SVMcon along with several other predictors during the CASP7 experiment.

Comparison with CMAPpro on the same Benchmark
SVMcon is trained to predict medium-to long-range contacts (sequence separation >= 6) as in [36], which are not captured by local secondary structure. We train SVMcon on the same dataset used to train CMAPpro [28,35] and test both programs on the same test dataset. The training dataset contains 485 proteins and the test dataset contains 48 proteins. The sequence identity between the training and testing datasets is below 25%, a common threshold for ab initio prediction [36].
We use sensitivity and specificity to evaluate the performance of SVMcon and CMAPpro. Sensitivity is the percentage of native contacts that are predicted to be contacts. Specificity is the percentage of predicted contacts that are present in the native structure. The contact threshold is set at 8 Å between Ca atoms. The sensitivity and specificity of a predictor depend also on the threshold used to separate 'contact' from 'non-contact' predictions. To compare SVMcon and CMAPpro fairly, we choose to evaluate them at their break-even point, where sensitivity is equal to specificity as in [37]. At the break-even point, the sensitivity and specificity of SVMcon is 27.1%, 4% higher than CMAPpro. Thus on the same benchmark dataset, SVMcon yields a sizable improvement.
We also compare the accuracy of SVMcon with the random uniform baseline algorithm consisting of random independent coin flips to decide whether each residue pair is in contact or not. Since the medium-and long-range contacts account for 2.8% of the total number of residue pairs with linear separation >= 6, the probability for the coin to produce a contact is set to 2.8%. As a result, the sensitivity and specificity of the random baseline algorithm is 2.8% at the break-even point. Thus SVMcon yields a nine-fold improvement over the random baseline.
Since the contact prediction accuracy varies significantly with individual proteins and their structure classes [29], we calculate the contact prediction specificity (or called accuracy) and sensitivity (or called coverage) for each test protein (Table 1). For each one, we select up to L (protein length) predicted contacts ranked by SVM scores because the total number of true contacts is approximately linear to the protein length [24]. The results show that in many cases (e.g. 1QJPA, 1DZOA, 1MAIA, 1LSRA, 1F4PA, 1MSCA, 1IG5A, 1ELRA, 1J75A), the prediction accuracy and coverage are > 30%.
However, for some proteins such as 1SKNP, the prediction accuracy is pretty low. We observe that the contact prediction accuracy is related to the the quality of multiple sequence alignment, the prediction accuracy of secondary structure, and the proportion of beta-sheets. Consistent with previous research [29,37], the contacts within betasheets in beta, alpha+beta, and alpha/beta proteins are predicted with higher accuracy than the contacts between an alpha helix and a beta strand or between alpha helices. We think that the strong restraints between beta-strands such as hydrogen-bonding probably contribute to the improved accuracy. Figures 1 and 2 show the native 3D structure and the predicted contact map of a good example (protein 1DZOA), respectively. In this case, 2L (240) predicted contacts with sequence separation >= 6 are selected. The contact prediction accuracy and coverage are 39.2% and 42.5%, respectively. It is shown that the predicted contact clusters ( Figure 2) recall most native beta-sheet pairing patterns of the protein (Figure 1). And it is interesting to see most false positive contacts are also clustered around the true contacts. Thus, these noise may not be very harmful during the process of reconstructing protein structure from the contacts. Furthermore, to investigate the relationship between the SVM contact map predictions and the structure classes, we compute the average accuracy and coverage of contact predictions in the six SCOP [38] structure classes ( Table 2). The contact prediction accuracy of proteins having betasheets (alpha+beta, alpha/beta, beta) is higher than that of alpha helical proteins, which is consistent with previous observations [29]. According to Table 2, the average coverage is about 20% and the accuracy ranges from 21 to 37%. This level of accuracy is probably good enough (or at least helpful) for constructing an ab initio low-resolution structure, since previous experiments show that only L/5 or even less residues contacts are required to reconstruct a low resolution structure for a small protein [5,8,[39][40][41][42], taking into account the inherent physical restraints of protein structures. However, the challenge is to develop algorithms to reconstruct a protein structure from a noisy predicted contact map, where contact restraints are much less reliable than the experimental contacts determined by NMR techniques.

Comparison with seven other Predictors during CASP7
SVMcon participated in the CASP7 experiment in 2006 and was evaluated along with seven other contact map Column 1 lists the protein name (PDB code + chain id). The chain id of a single-chain protein is denoted by "A" instead of "-". Column 2 lists chain lengths, ranging from 46 to 198. Column 3 lists the SCOP structure class, alpha, beta, a+b, a/b, small, and coil-coil represent six SCOP protein classess (all alpha helix, all beta sheet, alpha helix + beta sheet, alpha helix alternating beta sheet, small protein, and coil-coil), respectively. Columns 4 and 5 report the prediction accuracy (specificity) and coverage (sensitivity) for each protein. Accuracy is the number of correct predictions divided by the total number of predictions. Coverage is the number of correct predictions divided by the total number of true contacts. The raw number of correct preditions, all predictions, and true contacts are also reported in the brackets.
predictors. The CASP7 evaluation procedure focuses on inter-residue contact predictions with linear sequence separation >= 12 and >= 24 respectively [15]. Up to L/5 of the top predicted contacts were assessed, where L is the length of the target protein. These evaluation metrics are also similar to those used in the past editions of the Critical Assessment of Fully Automated Structure Prediction Methods [43][44][45] and in the EVA contact evaluation server [46]. We use the similar procedure to compute accuracy (specificity) and coverage (sensitivity) for all server predictors.
We compare SVMcon with the other contact map predictors on the 13 out of 15 CASP7 de novo domains whose structures have been released. The contact map predictors participating in CASP7 include SVMcon, BETApro [37], SAM-T06 [47], PROFcon [32], GajdaPairings, Distill [34,48], Possum [19], and GPCPRED [29]. Their contact predictions were downloaded from the CASP7 website. Overall, these results on the CASP7 dataset show that the accuracy and coverage of protein contact prediction are still low. However, these results are an important step towards reaching the milestone of an accuracy level of about 30%, required for deriving moderately accurate (low resolution) 3D protein structures from scratch [5,8,[39][40][41][42] (Also, Dr. Yang Zhang, personal communication at the CASP7 conference). In particular, these predictors tend to predict different correct contacts. Thus, a consensus combination of contact map predictors may yield more accurate contact map predictions, which in turn could significantly improve 3D structure reconstruction. Since the stakes of contact map prediction are high, a community-wide effort for improving contact map prediction should be worthwhile (Dr. Burkhard Rost's lecture slides at Columbia University).
It is also worth pointing out that the CASP7 de novo dataset is too small to reliably estimate the accuracy of the predictors. So one should not over-interpret these results. Indeed, when we use a larger CASP de novo dataset of 24 domains classified by Dr. Dylan Chivian from Dr. David Baker's group to evaluate the predictors (results not shown), the accuracy of SVMcon and BETApro are very close for both sequence separations >= 12 and 24, and both remain among the top predictors.

Conclusion
We have described a new contact map predictor (SVMcon) that uses support vector machines to integrate a large number of useful information including profiles, second-  ary structure, solvent accessibility, contact potentials, residue types, segment window information [24,32], and protein-level information [32]. The method yields a 4% improvement over the state-of-the art contact map predictor CMAPpro. In the blind CASP7 experiment, SVMcon is ranked as one of the top contact predictors. The method represents an effort toward a good 2D structure prediction. It can be used to improve ab initio structure prediction [4] and analogous fold recognition [2,3]. The web server, software, and source code are available at the SVMcon website [49].

Data Sets
In the comparison with CMAPpro [28,35], we use the same training and testing datasets. The datasets are redundancy reduced. The pairwise sequence identity of any two sequences is less than 25%. The training and testing datasets contain 485 sequences and 48 sequences respectively.
We use PSI-BLAST to search each sequence against the NCBI non-redundant database and generate a multiple sequence alignment. The number of PSI-BLAST iterations is set to 3. The e-value for selecting a sequence is set to 0.001. The e-value for including a sequence into the construction of a profile is set to 10 -10 . Multiple sequence alignments are converted into profiles, where each position is associated with a vector denoting the probability of each residue type.
We extract only medium-and long-range residue pairs with sequence separation >= 6 as in [32], which are not captured by local secondary structures. We use a 8 Å distance threshold between Ca atoms to classify residue pairs into two categories: contact (positive, < 8 Å) or non-contact (

Input Features
We extract five categories of features for each pair of residues at positions i and j to evaluate their contact likelihood. In addition to the new features (e.g. pairwise information features), the choice of most features combines ideas from our previous contact map predictors, disulfide bond predictors [33,50], and beta sheet topology predictors [37], and from the PROFcon [32], the best predictor in CASP6.

Local window feature
We extract local features using a 9-residue window centered at each residue in each residue pair. For each position in the window, we use 21 inputs for the probabilities of the 20 amino acids plus gap, computed from multiple sequence alignments, 3 binary inputs for secondary structure (helix: 100, sheet: 010, coil: 001), 2 binary inputs for relative solvent accessibility (exposed: 10, buried: 01) at 25% threshold, one input for the entropy (-∑ k p k logp k ) as a measure of local conservation. Here p k is the probability of occurrence of the k-th residue (or gap) at the position under consideration. Secondary structure and relative solvent accessibility are predicted using the SSpro and ACCpro programs in the SCRATCH suite [27,35,51]. Thus the two local windows produce 2 × 9 × 27 features.

Pairwise information features
For each pair of positions (i, j) in a multiple sequence alignment, we compute the following features. One input corresponds to the mutual information of the profiles of The key contacts within anti-parallel strand pairs (3,4), (4,5), and (5,6) are recalled. A few contacts within the parallel strand pair (1,2) are also predicted correctly. However, very long range contacts between alpha helices and beta sheets are not predicted. And there are some false positives. It is interesting to see that most false positives are close to the true contacts. Thus, they may not be very harmful when being used as distance restraints to reconstruct protein 3D structure. Thus some information about correlated mutations is used in the inputs. We also use three inputs to represent Levitt's contact potential [52], Jernigan's pairwise potential [53], and Braun's pairwise potential [54] for the residue pairs in the target sequence.
Central segment window feature Central segment window corresponding to a window centered at position Q(i + j)/2)N has been shown to be useful for predicting whether the residues at position i and j are in contact or not [24,32]. We use a central segment window of size 5. For each position in the window, we use the same 27 features as the local window features. So the total number of features for the central window is 5 × 27. We also compute the amino acid composition (21 features), secondary structure composition (3 features), relative solvent accessibility composition (2 features) in the central segment window. The sequence separation (|i -j + 1|) for residue pair (i, j) are classified into one of 16 length bins (< 6, 6, 7, 8, 9, 10, 11, 12, 13, 14, < 19, < 24, <= 29, <= 39, <= 49, >= 50) using a binary vector of length 16, as in [32].
x y x y Column 1 lists six structure classes. Column 2 lists the number of proteins in each class. Other columns reports the accuracy and coverage of contact predictions in each category. The statistics is computed for sequence separation >= 6, 12, and 24, respectively. The last row reports the average performance on all 48 proteins. The accuracy of a+b and a/b is slightly higher than that of beta proteins, which is in turn higher than that of alpha proteins. The performance on small proteins (mostly alpha helical) lies between proteins containing beta-sheets (a+b, a/b, and beta) and alpha helical proteins. There is only one coil-coil protein, which does not have native contacts with sequence separation >= 24. The eight contact map predictors are evaluated on the 13 de novo domains of CASP7. The 13 domains include (T0296, T0300, T0307, T0309,  T0314, T0316 domain 2, T0319, T0347 domain 2, T0350, T0353, T0361, T0356

Protein information features
We also compute the global amino acid composition (21 features), secondary structure composition (3 features), and relative solvent accessibility composition (2 features) of the target sequence. In addition, we classify sequence lengths into four bins (<= 50, <= 100, <= 150, and > 150) using a binary vector of length 4, as in [32].
The detailed methods of generating features are described in the additional files [see Additional file 1, 2, 3].

Feature Selection
Feature selection is useful to improve the performance of machine learning methods, particularly when there is a large number of features as in this study. However, since there are more than 310,000 training data points, it takes about 12 days to conduct a round of training and testing on a Pentium-IV computer. Thus a thorough feature selection is currently not feasible. So we tried only to remove some features (pairwise profile correlation, pairwise mutual information, residue type, and protein information features) once a time to test how they affect the prediction accuracy. We find that removing these features slightly improve the accuracy by about 0.2%. However, it is not clear if the improvement is due to the random variation or due to the removal of the features. But at least, these features are not essential or being compensated by other similar features. Thus, a more thorough feature selection should be conducted to improve the performance when more computing power is available.

SVM Learning
For an input feature vector associated with a pair of residues, we use Support Vector Machines (SVMs) to predict if the two residues are in contact (positive) or not (negative). SVMs provide a non-linear classifier model by nonlinearly mapping the input vectors into a feature space and using linear methods for classification in the feature space [55][56][57][58]. Thus SVMs, and more generally kernel methods, attempt to combine the advantages of both linear and nonlinear methods by first embedding the data into a feature space equipped with a dot product and then using linear methods in the feature space to perform classification or regression tasks based on the Gram matrix of dot products between data points. A key property of kernel methods is that the embedding does not need to be given in explicit form, the Gram matrix of dot products K (x, y) = φ (x)·φ (y) between data points is all is needed to proceed with classification or regression. Here x and y are input data points, φ is the mapping from input space to feature space, and K is the kernel or similarity measure between the original data points. Given a set of training data points S = S + ∪ S -, where S + (resp. S -) represent the positive (resp. negative) examples, using the theory of structural risk minimization [55][56][57][58] We use SVM-light [59][60][61] to implement SVM classification on our data. We experimented with several common kernels including linear kernels, Gaussian radial basis kernels (RBF), polynomial kernels, and sigmoidal kernels. In our experience, on this data the RBF kernel K (x, y) = (or ) works the best. Using the RBF kernel, f (x) is actually a weighted sum of Gaussians centered on the support vectors. Almost any separating boundary or regression function can be obtained with such a kernel [62], thus it is important to tune the SVM parameters carefully in order to achieve good generalization performance and avoid overfitting.
We only adjust the width parameter γ of the RBF kernel, leaving all other parameters to their default value. γ is the inverse of the variance (σ 2 ) of the RBF and controls the width of the Gaussian functions centered on the support vectors. The bigger is γ, the more peaked are the Gaussians, and the more complex are the resulting decision boundaries [62]. After experimenting with several values of γ, we selected γ = 0.025.