DEEPre: sequence-based enzyme EC number prediction by deep learning

2.1 EzyPred

EzyPred (Shen and Chou, 2007) is a three-level EC number predictor, which predicts whether an input protein sequence is an enzyme, and its main class and subclass if it is. It uses two features, pseudo PSSM (Pse-PSSM) and FunD encoding. Pse-PSSM is developed from the pseudo amino acid, a highly innovative manually designed feature which has already been proved successful in a number of problems (Chou, 2009; Hayat and Khan, 2012). It encodes the PSSM of proteins with different lengths using a uniform length matrix, which not only preserves the average score of the amino acid residues in the whole sequence that were changed to a certain type of amino acid during the evolution process but also avoids the complete loss of the sequence order information. FunD encoding captures the local FunD information, which could be referred to Section 3.2.5. With these two features, EzyPred uses optimized evidence-theoretic k-nearest neighbor (OET-KNN) as the classifier, which is an improved version of KNN. By considering not only the label of the KNN of the input query data point but also the distance of the neighbors to the query data as the supporting evidence, OET-KNN alleviates the problem of the original version of KNN for being too sensitive to noise. Although having been developed for 10 years, EzyPred still remains as one of the state-of-the-art methods in predicting enzyme function. Its server is easy-to-use with a user-friendly interface as well.

2.2 SVM-prot

SVM-Prot was proposed in 2004 (Cai et al., 2003, 2004) and updated in 2016 (Li et al., 2016). It can not only predict enzyme functional families but also non-enzyme functional families. It represents the protein sequence using 13 properties, including AAC, polarity, hydrophobicity, surface tension, charge, normalized Van der Waals volume, polarizability, secondary structure, solvent accessibility, molecular weight, solubility, number of hydrogen bond donors in side chain and number of hydrogen bond acceptors in side chain. Employing composition, translation and distribution to encode each of the above properties, SVM-Prot can make prediction irrespective of sequence similarity. Specifically, composition specifies the fraction of amino acids with a particular property; translation specifies the transition percentage of one amino acid with particular property to another amino acid with different properties; distribution specifies the distribution of amino acids with certain property within the first 25, 50, 75, and 100% of the sequence. The original version used support vector machines (SVM) as the classifier, while the 2016 update made two more classifiers, KNN and probabilistic neural networks, available.

2.3 COFACTOR

COFACTOR (Roy et al., 2012; Zhang et al., 2017) is a structure-based protein function annotation web-server. In terms of EC number prediction, for an input structural model, which can be obtained either by experiments or computational modeling, it threads the structure against the template library, whose entries’ annotation has already been validated by experiments, to identify the template enzyme with the most similar folds and functional sites. Obtaining the template and assuming that structures determine function, the server assigns the EC number of the template enzyme to the query, with the confidence being evaluated by a function considering both the global similarity and the local similarity. In addition to enzyme function prediction, the server can predict the Gene Ontology (GO) terms and protein-ligand binding interactions as well. COFACTOR has been proved successful in protein–ligand binding interaction prediction in the CASP9 competition (Moult et al., 2011).

2.4 EFICAz

EFICAz (Arakaki et al., 2009; Kumar and Skolnick, 2012; Tian et al., 2004) is an EC number prediction server using combined approaches. In addition to using the sequence similarity, it also incorporates the PROSITE and PFAM database information. The original version consists of four components: (i) pairwise sequence comparison-based enzyme function inference, (ii) conservation controlled hidden Markov model (HMM) iterative procedure for enzyme family classification-based functionally discriminating residue identification, (iii) multiple PFAM-based functionally discriminating residue recognition and (iv) multiple PROSITE pattern recognition. Those four components work independently, determining the final prediction by voting. In the later updates in 2009 and 2012, two more components, multiple PFAM family-based SVM evaluation and conservation controlled HMM iterative procedure for enzyme family classification-based SVM evaluation, and larger databases were added. Although it is unable to make EC number annotation if the query sequence has no homology, this server works pretty well in practice with completely four digits assigned.

2.5 Deep learning and hierarchical classification

Since (Krizhevsky et al., 2012), deep learning has become an extremely popular machine learning method. Its two main architectures, convolutional neural network (CNN) and recurrent neural network (RNN), have made a profound contribution to many bioinformatic problems, such as genetic analysis (Xiong et al., 2015), sequence binding specificity prediction (Alipanahi et al., 2015), and cryo-EM image processing (Wang et al., 2016a). Instead of being a pure classifier that depends on the manually designed features such as SVM, CNN is considered as an end-to-end wrapper classifier, being able to perform feature extraction based on the classification result and improve the performance in a virtuous circle. As a complement to CNN’s capability of capturing significant features from a 2D or 3D matrix, RNN has the potential of encoding long term interactions within the input sequence, which is usually a 1D vector, such as the encoding of English words. In our article, we combined the advantages of CNN and RNN, using CNN to conduct feature extraction and dimensionality compression starting from the raw 2D encoding matrices, and using RNN to extract the sequential, long-term interactions within the input sequence.

A classification problem with a tree structure in the label space, such as the enzyme function prediction problem discussed in this article, is often regarded as a hierarchical classification problem. Because this kind of problems can be regarded as multi-label classification and multi-class classification at the same time, the solutions to the problem can be classified into three categories based on different angles to the problem (Silla and Freitas, 2011), namely, flat classification approach, local classifier approach, and global classifier approach. According to the property of our problem, we chose the local classifier approach, which constructs one classifier for each internal node, to be the overall strategy.

3.1 Datasets

We adopt three datasets in this paper. The first dataset is a widely used one from (Shen and Chou, 2007), constructed from the ENZYME database (released on May 1, 2007), with 40% sequence similarity cutoff. More details of that dataset could be referred to (Shen and Chou, 2007). This dataset is denoted as the KNN dataset in the rest of the paper.

Following the same procedure of constructing the KNN dataset, we constructed a larger dataset using up-to-date databases. The steps of constructing the dataset are as follows:

1. The SWISS-PROT (released on September 7, 2016) database was separated into enzymes and non-enzymes based on the annotation.

2. To guarantee the uniqueness and correctness, enzyme sequences with more than one set of EC numbers or incomplete EC number annotation were excluded.

3. To avoid fragment data, enzyme sequences annotated with ‘fragment’ or with <50 amino acids were excluded. Enzyme sequences with more than 5000 amino acids were also excluded.

4. To remove redundancy bias, we used CD-HIT (Fu et al., 2012) with 40% similarity threshold to sift upon the raw dataset, resulting in 22 168 low-homology enzyme sequences.

5. To construct the non-enzyme part, 22 168 non-enzyme protein sequences were randomly collected from the SWISS-PROT (released on September 7, 2016) non-enzyme part, which were also subject to the (ii–iv) steps.

This larger dataset would be referred to as the NEW dataset in the rest of this article.

Other than KNN and NEW, which will be used as the benchmark to evaluate the proposed method based on cross-fold validation, it is also important to test the generalization power of the proposed method. This can be done by training the model on one dataset, and testing it on an independent and non-overlapping dataset, to avoid being overfitted on a particular dataset. Thus, the third dataset, the benchmark dataset from (Roy et al., 2012), is used for cross-dataset validation. This non-homologous dataset was collected from PDB, satisfying two requirements: (i) the pair-wise sequence similarity within the dataset is below 30%, and (ii) there is no self-BLAST hit within the dataset to ensure that there are no enzymes that are homologous to each other in this set (Roy et al., 2012). All enzymes in this dataset have experimentally determined 3D structures. To avoid overlaps between the training and testing datasets, sequences contained in both our training dataset and this dataset were removed, which reduced the size of the dataset from 318 to 284. This benchmark dataset would be referred to as the COFACTOR dataset in the following. Table 1 summarizes the three datasets.

3.2 Sequence representation

The deep learning framework explained in Section 3.3 eliminates the necessity of performing manual dimensionality uniformization and building complex, manually designed features, which are unlikely to sustain the increasing amount and complexity of data, by conducting feature reconstruction and classifier training simultaneously. Therefore, we use the following raw features, constructed from the input sequence directly, to represent the sequences. Based on their dimensionality, they can be classified into two categories, sequence-length-dependent features and sequence-length-independent features. The first four features described below belong to the former while the last one belongs to the latter.

3.3 Classification model

The enzyme function prediction problem has a tree-structured label space, which makes it a typical hierarchical classification problem. To solve this kind of problems, we propose a level-by-level prediction framework, building a model for each internal label node. The model contains two main components, namely, the problem-specific feature extractor, which is able to perform dimensionality uniformity and feature extraction, and the classifier. Such a novel, end-to-end model can perform feature selection and classifier training simultaneously in a virtuous circle, making it more likely to achieve high performance.

4.1 Evaluation criteria

For the enzyme or non-enzyme prediction, since it is a binary classification problem, we use accuracy, Cohne’s Kappa Score (Viera and Garrett, 2005), precision, recall and F₁ score to evaluate the classifiers’ performance. For other predictions, since they are multi-class classification problems, we use accuracy, Cohen’s Kappa Score, Macro-precision, Macro-recall and Macro-F₁ score to evaluate the classifiers’ performance, whose definitions are in Supplementary Section S3.

4.2 Compared methods

For the cross-fold validation, in which training and testing are based on different parts within the same dataset, we compare our method with five other methods, including two state-of-the-art methods, EzyPred (Shen and Chou, 2007) and SVM-Prot (Li et al., 2016), and three baseline methods. One of the baseline methods uses SVM with the raw features used in our model; another baseline method uses SVM with Pse-PSSM; and the last baseline method uses the traditional neural network with our raw features. Due to the unchangeable database of EFICAz (Kumar and Skolnick, 2012) and COFACTOR (Zhang et al., 2017), we do not include them in the cross-fold validation comparison. However, we perform cross-dataset validation, where the training and testing are performed on different datasets, to compare our method with EzyPred, SVM-Prot, COFACTOR and EFICAz.

4.3 Cross-fold validation

Here we report the 5-fold cross validation results, which are shown in Figure 2. Our method almost always outperforms the other methods in both the KNN dataset and the NEW dataset across the five criteria and across the three hierarchical levels of prediction. As for the NEW dataset, DEEPre outperforms the other five methods consistently in Levels 0 and 1 prediction across the five criteria. As for the Level 2 prediction, the only criterion that DEEPre does not improve over the existing methods is the Macro-Precision, which is an unweighted average of precision of each label. The appearance of very small classes (e.g. subclass 1.20 only has 10 enzymes) in the second level prediction might be the reason for this result. In terms of the KNN dataset, although the smaller dataset makes the improvement of DEEPre over the other methods in Level 0 prediction less significant, it still significantly outperforms the other methods in Levels 1 and 2 classification.

4.4 Feature importance analysis

It is believed that both global features and local features determine the function of a protein. For detailed function, local information would weigh even more in determining it. The features extracted by the convolutional component and the recurrent component from PSSM and sequence raw encoding could be considered as global features while the FunD encoding would be considered as a local feature. We remove the three input raw encoding one by one and show the comparison of their performance on the NEW dataset. The comparison is shown in Figure 3A. It is clear that as the level goes deeper, the importance of FunD is evidently increasing, which demonstrates the well-recognized hypothesis. To further prove it, we design another experiment, in which we input more local feature encodings, including secondary structure and solvent accessibility, into our model. Details of this experiment could be referred to Supplementary Section S4. Figure 3B shows the performance comparison of this model and the previous model in Level 2 prediction. It is clear that the additional local features further improve the performance of our model, with accuracy improved by 1.8% while Macro-precision, Macro-recall and Macro-F1 score improved by at least 11%.

4.5 Cross-dataset validation

In this experiment, we directly compare the performance of different servers in predicting the first digit and the second digit of an enzyme. We use the COFACTOR benchmark dataset, which is proved to be a difficult dataset in the enzyme function prediction field (Roy et al., 2012), as the test dataset. First, we eliminate the sequences in the COFACTOR benchmark data which overlap with the DEEPre’s training database (NEW) by 40% sequence similarity filtering, reducing the data size from 318 to 284, to ensure that there is no bias in the DEEPre’s results. Then we input the remaining sequences to each server manually and collect the prediction results. For COFACTOR, since it is quite time-consuming to run the server, about 4 h to obtain the result for one query, we report the results from the original paper. As shown in Figure 4, for the first-digit prediction, DEEPre outperforms the other servers consistently across the five criteria, improving the accuracy by at least 6% over the other servers, including COFACTOR. This is significant because COFACTOR requires 3D structures of enzymes whereas DEEPre only requires the sequence information. On the other hand, we should admit that we have changed the original COFACTOR dataset to some extent by reducing the overlap between the training and testing sets, which might explain some of the performance difference between COFACTOR and DEEPre. We should also notice that all of those five servers have different training datasets, but those training datasets highly overlap with each other and each method was optimized on its corresponding dataset. In addition, although we removed overlapping enzymes from our training set and the COFACTOR test set, there may still be homologs of the enzymes in COFACTOR in our training set. However, DEEPre is a sequence-based statistical method, which explores the statistical properties of training data and does not benefit from knowing enzyme structures. Therefore, although the performance of DEEPre on COFACTOR may still be biased, the influence is not as much as that by nearest neighbor-based methods. We also perform comparison of different servers’ performance on the second-digit prediction (Supplementary Section S5). The results show that DEEPre and EFICAz both perform well on the second-digit prediction on the COFACTOR dataset. It is worth noting that EC numbers have regular corrections, such as deletions and transfers. We check all the corrections that are related to the test enzymes in the COFACTOR dataset and find that none of them influences the comparison reported here.

4.6 Third-digit and fourth-digit prediction

Using the same framework described above, we are also able to predict the enzyme’s third digit, which represents its sub-subclass, on the NEW dataset. The accuracy across all the sub-subclasses is 0.9415; the Kappa score is 0.8918; the macro-precision is 0.8942; the macro-recall is 0.8578; and the macro-F1 score is 0.8665. Regarding the fourth-digit prediction, more data are needed to perform normal machine learning training-and-testing procedure. For example, within the sub-subclass 1.1.1 in the NEW dataset, there are 188 classes. Each of those classes has <40 enzyme sequences, with 175 classes having <10 enzyme sequences. Using the current dataset with such distribution would lead to unreliable results.

4.7 Case study

Glutaminase is a phosphate-activated enzyme, which catalyzes the first step of glutaminolysis, hydrolysing glutamine into glutamate (Curthoys and Watford, 1995). The alternative splicing of its messager RNA results in its three isoforms, with Isoforms 1 and 3 being capable of catalyzing while Isoform 2 lacking the catalytic activity (Li et al., 2017). To validate our model’s ability to distinguish the different functionality of different isoforms, we obtained the sequences of the three Glutaminase isoforms from the UniProt and put them into our model. Our model predicted that Isoforms 1 and 3 of Glutaminase were hydrolases acting on carbon-nitrogen bonds, being consistent with the experimental results. Our model also recognized Isoform 2 as non-enzyme, which is consistent with the experimental result as well.

Aurora kinases B is a key enzyme regulating chromosomal segregation during mitosis, ensuring correct chromosome alignment and segregation as well as chromatin-induced microtubule stabilization and spindle assembly (Carmena et al., 2009). Over-expression of it is possible to cause unequal distribution of genetic information, resulting in aneuploid cells, which may become cancerous (Sorrentino et al., 2005). Aurora kinases B has five isoforms resulted from alternative splicing. Four of them have roughly equal length with high similarity, while Isoform 3, having high expression in the metastatic liver with no expression in the normal liver, is only half of the length of the” canonical” isoform (142 amino acids versus 344 amino acids). Despite its much shorter length, the isoform does not lose its functionality. To further validate our model’s ability of handling isoforms’ functionality prediction, we collected the sequence of the five isoforms from the database and put them into our model. Our model’s result is consistent with the experimental results. Particularly, our model predicted the functionality of the Isoform 3 successfully, despite its sequence’s large difference from the ‘canonical’ sequence.

The detailed performance comparison of different servers on these two case studies could be referred to Supplementary Section S6. Among the five compared methods, only our method and EzyPred produced correct predictions for both cases.