Accurate classification of membrane protein types based on sequence and evolutionary information using deep learning

Background Membrane proteins play an important role in the life activities of organisms. Knowing membrane protein types provides clues for understanding the structure and function of proteins. Though various computational methods for predicting membrane protein types have been developed, the results still do not meet the expectations of researchers. Results We propose two deep learning models to process sequence information and evolutionary information, respectively. Both models obtained better results than traditional machine learning models. Furthermore, to improve the performance of the sequence information model, we also provide a new vector representation method to replace the one-hot encoding, whose overall success rate improved by 3.81% and 6.55% on two datasets. Finally, a more effective model is obtained by fusing the above two models, whose overall success rate reached 95.68% and 92.98% on two datasets. Conclusion The final experimental results show that our method is more effective than existing methods for predicting membrane protein types, which can help laboratory researchers to identify the type of novel membrane proteins.

According to their functions, membrane proteins can be classified into three classes: integral, peripheral and lipidanchored [9]. Based on the direct interaction relationship between membrane proteins and lipid bilayers, the three classes can be further extended into eight basic types: (1) type I membrane proteins, (2) type II membrane proteins, (3) type III membrane proteins (4) type IV membrane proteins, (5) multipass transmembrane proteins, (6) lipid chain-anchored membrane proteins, (7) GPI-anchored membrane proteins, and (8) peripheral membrane proteins. Among them, Types I, II, III and IV are singlepass transmembrane proteins, and detailed descriptions of their differences are given in [9].
The most popular feature extraction methods for predicting membrane protein types are based on sequence information. The pseudo-amino acid composition (PseAAC) [17] method incorporates information about sequence order. Local amino acid composition (LAAC), local dipeptide composition (LDC), global descriptor (GD), Lempel-Ziv complexity (LZC), autocorrelation descriptor (AD), sequence-order descriptor (SD) and Hilbert-Huang transform (HHT) are proposed based on amino acid classification and physicochemical properties [18], and these methods are partially applied in our previous research [19]. The peptide composition method, such as dipeptide composition (DipC) [20] and tripeptide composition (TipC) [21] are also powerful sequence information-based feature extraction methods.
In this paper, we improve the performance in membrane protein type prediction. The main contributions of the paper are summarized as follows: First, we present two deep neural network (DNN) models to process the sequence information and evolutionary information separately. These models both achieve better performance than traditional machine learning models. Second, by using the DNN model with convolutional and recurrent layers, we can remove the burden of feature engineering and the reliance on domain experts. Third, we provide a new vector representation based on autoencoder and physicochemical property indexes. Numerous experimental results prove that the new vector has more powerful representation ability than one-hot encoding. Finally, a more effective model is constructed by fusing the above two models with the ensemble learning.
The DNN model has shown its superiority in the field of bioinformatics [28][29][30][31][32][33][34][35][36][37]. However, to the best of our knowledge, this paper is the first to propose DNN models for the prediction of membrane protein types. The two DNN models proposed in this paper can process the sequence information and evolutionary information, separately. When processing sequence information, we use a 1D convolutional layer [38] and bidirectional long short term-memory (Bi-LSTM) layers [39], and when processing evolutionary information, a 2D convolutional layer and CapsNet layer [40] are adopted. Among these layers, compared to convolutional network, CapsNet is a new network that can extract local position information.

Datasets
In this work, we use two datasets to evaluate the performance of our method, hereafter referred to as Dataset 1 and Dataset 2. Dataset 1 consists of two benchmark sets, the training set and the testing set, which are both from paper [41] and used in previous studies [9,10,12,13,20]. Dataset 2 is constructed by Chen [12]; it also consists of two benchmark sets. Table 1 shows the concrete distribution of Dataset 1 and Dataset 2.

Architecture of the proposed DNN
To improve the performance of the prediction of membrane protein types, we build the DNN model with a keras framework (http://www.keras.io). Considering the respective characteristics of sequence information and evolutionary information, we build two different models to process them, which are hereafter referred to as the sequence information model and evolutionary information model, respectively.

Sequence information model
In this model, the input sequences are converted into numerical vectors of fixed length by truncating the long sequence and padding the short sequence with '0' . The length of the sequence ranges from 50 to 5000 in our dataset. Although the length of the sequence varies Table 1 The distribution of samples for Dataset 1 and Dataset 2   Membrane protein types  Dataset 1  Dataset 2   Training set  Testing set  Training set  Testing set   Type I  610  444  561  245   Type II  312  78  316  79   Type III  24  6  32  9   Type IV  44  12  65  17   Mutipass  1316  3265  1119  greatly, sequences with a length of less than 1500 account for almost 98% of the total dataset. Therefore, taking the utilization of information and the computational complexity into account, we decided to set the fixed length to 1500. Each of the basic 20 amino acids is converted into a number of 0-19. Ultimately, to prevent different encoding from affecting the performance of model, one-hot encoding is used before inputting the vector into the model. The architecture of the sequence information model, as shown in Fig. 1, consists of one 1D convolutional layer (filters: 256, kernel__size: 15, strides: 10, padding: 'same' , activation: 'relu'), two Bi-LSTM layers (units: 128, dropout: 0.5 and rest default setting) and one fully connected layer (softmax). The 1D convolutional layer can reduce the complexity of the Bi-LSTM layer that reduces the shape of the output matrix from 1500×256 tensor to 150×256. Moreover, the 256 filters in the layer can also extract rich information from the sequence. Next, two Bi-LSTM layers are applied, which can identify features separated by large gaps. In previous studies, a normal LSTM layer is mainly used to process sequence information [29,30]. Compared to normal LSTM, the Bi-LSTM Fig. 1 The sequence information model uses 1D Conv and Bi-LSTM layers. First, the membrane protein sequences are encoded into a matrix of size 1500×20. Then, these matrixes are fed to the 1D Conv layer (filters: 256, kernel size: 15, strides: 10, padding: same, activation: relu). Next, two Bi-LSTM layers (units: 128, dropout: 0.5 and the remainder as default settings) are used to identify features separated by large gaps. Last, the output from the Bi-LSTM layer is passed through a final dense layer that uses a softmax function to obtain the probability value belonging to each membrane protein type layer we used can take advantage of both historical input and future information, which performs better for our task. Each Bi-LSTM unit comprises three gates, which are shown in Fig. 1. Their calculation process is as follows: hidden cell states (6) where W f , W i , W C and W o denote the weight matrix of the forget gate, input gate, candidate cell states and output gate, respectively, and b f , b i , b C , and b o denote the bias of the forget gate, input gate, candidate cell states and output gate, respectively. The forget gate decides which information is discarded from the cell state, where h t−1 represents the output of the previous cell, and x t represents the input of the current cell. σ represents the sigmod function. The input gate decides how much new information is added to the cell states, and it uses two steps to accomplish this task: First, a sigmoid layer determines which information needs to be updated; a tanh layer generates a vectorC t , which is the alternative content to update. Then, we combine two parts to update the state of the cell. The output gate will determine what value to output. First, we run a sigmoid layer to determine which part of the cell state will be output. Next, we process the cell state through tanh (obtaining a value between -1 and 1) and multiply it by the output of the sigmoid gate. Eventually, we will only output the part of our output that we determined. Then, the output from Bi-LSTM layer is passed through a final dense layer that used a softmax function to obtain the probability value belonging to each membrane protein type. We compile our keras model with the ' Adam' optimizer (lr: 0.001, beta__1: 0.9, beta__2: 0.999, decay=0). Finally, we find that adding the BatchNormalization layers (default setting) after three mainly layers in our model helps to accelerate convergence.

Evolutionary information model
In this paper, evolutionary information mainly refers to the PSSM, which was first proposed by Jones for predicting the secondary structure of proteins [42]. The PSSM can discover protein sequences that have evolved relationships with search sequences. It is expressed as follows: where L is the length of the sequence, and M i→j denotes the score of the amino acid in the i-th position of the amino acid being mutated to j-th position of the amino acid during the evolution process. Because membrane protein sequences vary in size, their PSSM shapes are also different. Similar to the sequence information model, we also truncate the long sequence and pad the short sequence with '0' to convert the PSSM matrix into a 1500×20 matrix. The architecture of evolutionary information model is shown in Fig. 2. It consists of two 2D convolutional layers, two average pooling layers, one PrimaryCaps and one fully connected layer. The first 2D convolutional layer (filters: 256, kernel__size: (5,5), strides: 1, padding: 'valid' , activation: 'relu') is used to extract rich and effective features from the PSSM. Next, the first 2D average pooling layer (pool__size: (2,2) and the remainder as default settings) is used to retain the main features while reducing the parameters and calculations of the next layer to prevent overfitting. Although paper [40] argued that the pooling layer may throw away information about the precise position of the entity within the region, we find that the pooling layer and CapsNet work together to obtain the best results in our task. One possible reason is that the pooling layer can improve the robustness of extracted features. Then, the second 2D convolutional layers (filters: 128, kernel__size: (5,5), strides: 1, padding: 'valid' , activation: 'relu') and the second 2D average pooling layer (pool__size: (2,2) and the remainder as default settings) are reused to obtain more advanced features. Moreover, each convolutional layer has a dropout technique with the rate 0.5 to prevent the model from overfitting. The first four layers are designed to increase the representation power of CapsNet.
The PrimaryCaps layer contains 4 primary capsules that accept the basic features detected by the first four layers and generate a combination of features. The 4 main capsules of this layer are essentially similar to the convolutional layer. Among them, the 1×1 convolution kernel used in our study not only achieves cross-channel interaction and information integration, but also reduces the number of convolution kernel channels. Each capsule applies sixteen 1×1×64 convolution kernels to the 372×2×128 input tensor, thus generating a 372×2×64 output tensor. With a total of 4 capsules, the output is a 372×2 × 4×16 tensor. Since the length of a capsule Fig. 2 The evolutionary information model uses 2D Conv and CapsNet layers. First, the PSSM matrix is processed as the matrix with the same shape. These matrixes are fed to the 2D Conv layer (filters: 256, kernel size: (5,5), strides: 1, padding: valid, activation: relu) and average pooling layer (pool size: (2,2) and the remainder as default settings) to extract the rich and effective features. Next, the 2D Conv layer(filters: 128, kernel__size: (5,5), strides: 1, padding: 'valid', activation: 'relu') and average pooling layer(pool__size: (2,2)) are reused to extract more advanced features. The PrimaryCaps layer is the convolutional capsule layer, which has size 1×1 convolution kernels and 4 channels of 16D capsules. The ProteinCaps layer has eight 16D capsules to represent one of membrane protein types. Finally, the L2-norm of each capsule vector is calculated to indicate the probability of each type represents the probability that the entity presented, the CapsNet layer needs a new activation function, which is called the squashing function. It is a novel nonlinear activation function that accepts an input vector and then compresses its length to [0,1] as output, which is calculated as follows: where V j and s j are the vector output and input of the CapsNet, respectively. The last layer contains eight membrane protein type capsules, one for each type of membrane protein. The calculation between the PrimaryCaps layer and ProteinCaps layer is shown in Fig. 3. Among them, u i represents the ith capsule in PrimaryCaps. W ij is the weight matrix that represents the spatial relationship and other important relationships between low-level capsules and high-level capsules. There are eight capsules (V j , j ∈[ 1, 2, · · · , 8]) in ProteinCaps, each of which receives the inputs from all capsule outputs in PrimaryCaps. V j is calculated by the weighted sum ofũ ji and then is passed through the squashing function. Here, the weights c ij are obtained by an iterative dynamic routing process that can be found in [40].
The output of the ProteinCaps layer has eight 16D vectors. During training, for each training sample, the loss value is calculated as follows: where the value of T c is decided by the correct label. If the correct label corresponds to the membrane protein type in ProteinCaps, then T c = 1; otherwise, T c = 0 . For the value of the other hyperparameters m + , m − and λ, we used the default values 0.9,0.1 and 0.5, respectively, which are suggested in [40].Then, eight loss values are added to obtain the final loss.

Fusion information model
Many studies have proved that the fusing information method has better results than the single information method [12,20,[43][44][45]. Therefore, we also present a fusion information model based on ensemble learning. Among these methods, the ensemble learning method Computation between the PrimaryCaps and ProteinCaps layers.There are 2976 16D capsules (each u i is an 16D vector) in PrimaryCaps. Eacĥ u j|i is produced by multiplying u i by a weight matrix W i,j (16×16). Capsule V j (16D vector and j ∈ [1,8]) in ProteinCaps is produced by a weighted sum over allû j|i and the squashing nonlinear activation function. The parameter c i,j is determined by the iterative dynamic routing process that we adopted is a stacking algorithm, and its details can be found in our previous work [19]. Figure 4 shows the complete flow chart of the stacking algorithm, where the meta-classifier that we used is multiple logistic regression (MLR).
Research has shown that using the class probability output of the base classifier as the input of the meta classifier performs better [46]. Then, we use P(x) as the output of the base classifier, it can be represented as follows:

Model training
In most experiments, the DNN models are trained using identical training strategies. The dataset is divided into three parts, the training set, validation set and testing set, which play different roles in our work. The training set is used to train the model, and the validation set is used to adjust parameters and select the best model. The testing set is used to evaluate the performance of the model at the end. In our work, 20% of samples are separated from the original training set as a validation set. Then, the distributions of three sets in two datasets are shown in Table 2. All DNN models are implemented using keras 2.1.2. Model training and testing are performed on a personal computer equipped with an Nvidia GTX 1060 GPU.

Model evaluation
We use sensitivity (Se), specificity (Sp), accuracy (ACC) and overall success rate (OSR) to evaluate the classification performance. They are defined as follows: where TP, TN, FP and FN represent true positive, true negative, false positive and false negative, respectively. The Matthews correlation coefficient (MCC) is also used in our work. It is generally considered to be a more balanced indicator, even if the sample content of the two categories varies widely [47]. The MCC is essentially a correlation coefficient between the actual classification and the predicted classification. Its value range is [-1, 1]. A value of 1 indicates a perfect prediction for the subject. A value of 0 indicates that the predicted result is not as good as the result of random prediction, and -1 means that the predicted classification is completely inconsistent with the actual classification. The MCC is defined as follows: The construction of vector representation One-hot encoding ignores the possible relationship between amino acids. For example, alanine and cysteine have a certain relationship in their physicochemical properties, yet '0' and '1' do not provide much informtaion about this relationship. Moreover, we usually need more data to train because the amino acids are stored as sparse matrixes. Using a vector representation method can effectively solve the above problems. Continuous bag of words (CBOW) and Skim-Gram [48] are two popular vector representation methods that are widely used in the field of natural language processing (NLP). Limited by the size of protein samples and the characteristics of the protein sequence, these methods are no longer applicable. In this paper, we provide a new vector representation method based on autoencoder and the physicochemical properties of amino acids. This method can achieve a more compact representation of input symbols and yield semantically similar symbols close to each other in vector space. The AAindex is a database of numerical indexes representing various physicochemical properties of amino acids and pairs of amino acids [49]. The database now has a total of 557 indexes, which all come from the published paper. In this research, the physicochemical properties indexes with 'NA' are screened out, and the remaining 537 indexes are adopted. Using a physicochemical properties index that represents the amino acids is a good idea to replace the one-hot encoding. However, this method would consume a large amount of storage space and greatly increase the computational complexity. To process the above problems, autoencoder is used to create a more effective representation.
Autoencoder is a neural network that can capture the most important features of data. It consists of the input layer, multiple hidden layers and output layer. The input and output of the autoencoder are consistent, and the goal is to reconstruct itself using some high-level feature recombination. If the number of nodes of the hidden layer is smaller than the input layer and the output layer, it represents the same low-density information and is a centralized representation of the input data obtained from the learning. In our work, we take 537 physicochemical properties indexes of each amino acid as input, reconstruct this informatin with the autoencoder, and then obtain a more effective vector representation from the intermediate layer. The architecture of the autoencoder that we used is shown in Fig. 5.

The performance of vector representation
To explore the advantage of the new vector representation method, we compared the performance between the onehot encoding and our method in the sequence information model. We note that our process of parameter adjustment is all on the validation set and not reparameterized to apply on the testing set. Figure 6 shows the overall success rate of vector representations using different dimensions in Dataset 1 and Dataset 2. To verify whether the vector representation that we proposed is better than one-hot encoding, we also indicated the overall success rate of the one-hot encoding method in the figures. It can be seen from the figure that our method is better in most cases. Taking both the prediction performance and computational complexity into account, we finally chose the vector representation with 10-dimension to replace the one-hot encoding.  Furthermore, to verify whether the vector representation is helpful to improve the performance of the model, we also randomly generate a 10-dimensional vector in the interval [0,1] for comparison experiments; the results are shown in Table 3. From the table, we find that the vector representation that we proposed is more effective than other methods both on Dataset 1 and Dataset 2. In addition, we find that the random method is even worse than the one-hot encoding. Although the random method has a lower dimension, it cannot express the relationship between amino acids well, which makes the result worse. The final experimental results illustrate the reliability of the new vector representation method that we proposed.  Table 4 gives the overall success rate of the comparison between the sequence information model and traditional models on Dataset 1 and Dataset 2, where four popular features used for comparative experiments are AAC, DipC, TipC and PseAAC. Among them, AAC, DipC and Tipc are amino acid composition-based methods; PseAAC is the method based on amino acid composition and physicochemical properties. These features are all extracted from the membrane protein sequence. The machine learning algorithm used to compare random forests. Research has illustrated its effectiveness in the prediction of membrane protein types [50]. The optimal parameters in random forests are obtained by using the grid search method and 5-fold cross-validation on the training set. From the table, we find that our DNN model has a better overall success rate than all the machine learning methods. Compared with the DipC, which has the best results among the four sequence features, our model improved by 5.43% and 8.30% on Dataset 1 and Dataset 2, respectively. Tables 5 and 6 give the values of Se, Sp, ACC and MCC for eight membrane protein types using four traditional models and our DNN model on Dataset 1 and Dataset 2, respectively. Although our model is not as good as the traditional model in some indicators, our model achieves better results in most cases. Furthermore, we find that it     is difficult for the method combining feature engineering and the machine learning algorithm to predict the small types, especially in Type 4 membrane protein samples where the value of Sn and Mcc is zero in most traditional models. However, our model can predict the small type samples well. The reason may be the powerful feature extraction capabilities of the convolutional layer and Bi-LSTM layer that we used. The experimental results prove that our proposed sequence information model has better performance than traditional models.

The effective of evolutionary information model
Then, we use the same experimental method to compare the performance between the evolutionary information model and seven popular PSSM-based tradition methods proposed by other researchers, including EDPSSM, KPSSM, CPSSM, TriPSSM, RPSSM, CoPSSM and PsePSSM. They are all powerful features extracted from PSSM. Table 7 gives the overall success rate between seven traditional models and the evolutionary information model on Dataset1 and Dataset 2. From the table, we find the best method for the traditional model is TriPSSM, which utilizes PSSM linear probabilities to compute features. Although TriPSSM can obtain the outstanding grades on membrane protein types prediction, we find that it has 8000-dimensional features; the high dimensionality of features in turn increases the computing complexity. Our model is still effective since it can use GPU acceleration to yield better results in less time.  Tables 8 and 9 report the results, measured in Se, Sp, ACC and MCC over all methods on Dataset 1 and Dataset 2. From the tables, we can draw a conclusion that our model is more effective than the traditional model in evolution information extraction. The convolutional layer has a powerful ability to extract features from the evolutionary information, and the CapsNet layer can extract local position information. These abilities explain why our DNN model could achieve better results.

The effective of stacking ensemble method
To investigate on how the fusion model impacts the performance, we report the overall success rate of the validation set and testing set using two single models and    the fusion model (Table 10). From the table, we find that the fusion model achieves better results both on datasets, demonstrating its superiority in prediction. Figures 7 and 8 give the comparison results, measured in Se, Sp, ACC and MCC over three models on Dataset 1 and Dataset 2,respectively. We find that the fusion model achieves the best results in most cases. However, the performance of small-scale forecasts has not yet met our expectations. Other effective strategies for solving imbalance problems may contribute to improving our method.
Furthermore, the prediction of membrane protein types requires not only high classification accuracy but also certain operational efficiency. Table 11 lists the running time for different models when predicting 100 membrane protein samples. We find that the process of fusion takes approximately 1.93e-10s, and the running time of the fusion model is almost the sum of two single models. Though the time complexity of the fusion model is higher, it is also within the tolerance range.

Comparison with the existing methods
To demonstrate the performance of our DNN model in practical use, we compared our DNN model with eight state-of-the-art machine learning methods in the prediction of membrane protein types. These models all proceed through complex feature engineering, including feature extraction, dimension reduction and so on. Table 12 represents the overall success rate between our DNN model with eight models on Dataset 1 and Dataset 2. Taking the accuracy of prediction and the degree of domain experts dependence into account, we find the DNN model not only achieves the best performance both on Dataset 1 and Dataset 2, but also removes the burden of feature engineering, which illustrates the effectiveness of our model.

Conclusion
In this paper, we propose two DNN models to process the sequence information and evolutionary information separately. Among these models, 1D convolutional layers and Bi-LSTM layers are used for the sequence information model, while the evolutionary information model uses 2D convolutional layers and CapsNet layers. In the comparative experiments, we compare these models with the traditional model. Numerous experimental results demonstrate that our two proposed DNN models not only obtain more competitive results in membrane protein type prediction but also can remove the burden of feature engineering and the reliance on domain experts. In addition, to obtain a better prediction performance, the stacking ensemble method is used to fuse the above two models.
Furthermore, we provide a new vector representation method based on autoencoder and the physicochemical properties of amino acids to improve the performance of the sequence information model. Compared with the one-hot encoding, the new vector representation method not only represents the relationship between amino acids well but also can effectively improve the prediction performance.
Finally, our method can be applied not only to membrane protein type prediction but also to other fields of bioinformatics [51][52][53][54][55][56][57][58][59][60][61][62][63]. However, there is still room for further investigation. For example, the problem of an imbalanced dataset had a negative effect on the accuracy of small-sized types. Some strategies for processing with the imbalanced dataset may improve the performance of our DNN model.