MultiPep: a hierarchical deep learning approach for multi-label classification of peptide bioactivities

Defined bioactivity classes

Based on the downloaded peptides and their annotations, 20 bioactivity classes were defined (see Supplementary Information and Supplementary Table S1). The complete list of classes and how the used databases have contributed to their generation can be seen in Table 2. Some database class contributions were based on overlaps. For instance, the class “ACE inhibitor” stems only from the database BIOPEP-UWM, however, the database LAMP2 contributed with 39 peptides, because it had peptides from another class that overlapped with the peptides from the “ACE inhibitor” class.

Bioactivity classes were only created if 100 or more peptides were available for that class, and if the class had more than two peptides that were unique for that specific class. We believed that if the classes did not live up to these criteria, the classes would be to information-poor for a classifier to use. Some of the classes had major overlaps, while other classes had none (Supplementary Fig. S2). The largest class was the Antibacterial class with 14 362 peptides and the smallest was the Opioid class with 117 peptides. Analyses to see how the defined bioactivity classes differed in terms of overall peptide lengths and amino acids distributions were made (Supplementary Figs S3 and S4). Some of the bioactivity classes are very general (antimicrobial, peptide hormone, and antihypertensive) while other classes are more specific (ACE inhibitor, opioid, and dipeptidyl peptidase inhibitor). No negative class with peptides without bioactivity was defined, as the many different classes generated negative classes for each other. If required, the applied databases can be used to acquire more information about the classes.

Overall network architecture

MultiPep is a deep neural network implemented using Keras (http://www.keras.io) and Tensorflow [34] in Python. MultiPep’s architecture is inspired by a dendrogram from a hierarchical clustering of our defined bioactivity classes’ mutual overlaps (Fig. 1A and Supplementary Fig. S2). The max-normalized class overlap values were clustered using the complete linkage type and Euclidean distance as metric. We used the hierarchical clustering algorithm provided by SciPy [35]. The cutoff for the clustering was based on SciPy’s default cutoff value. Based on the classes’ divisions, a dendrogram template that consisted of five lower level class clades with two, six, three, five, and four classes in each, respectively, was created. Above the class clades, upper levels that connected the lower levels until all class clades were amalgamated were made (Fig. 1A).

Using the dendrogram template, our deep neural network was constructed (Fig. 1B). The network has “network class clades,” which consist of small convolutional neural networks (CNNs) and class-clade-specific output layers. See the section “Class-clade-specific convolutional neural networks” for a description of the CNNs. The network class-clade-associated output layers is connected to upper level output layers. Only the network class-clade-associated output layers (Output 4_1, Output 4_2, Output 4_3, Output 4_4, and Output 4_5 in Fig. 1B) contain weights. These layers are regulated by the sigmoid activation function. The upper level output layers (Output 1, Output 2_1, Output 2_2, and Output 3 in Fig. 1B) are concatenation layers that merge the two connected lower level output layers’ maximum output values. Table 3 shows which bioactivity classes the upper level output layers represent. To provide an example, if Output level_2_2_2 outputs 0.8 given a peptide, this peptides is predicted to belong to one of the bioactivity classes “antifreeze,” “cytokines/growth factors,” “antioxidative,” “drugdelivery,” “opioid,” “ACE inhibitor,” “antihypertensive,” “antidiabetes,” or “dipeptidyl peptidase inhibitor” by a score of 0.8.

Class-clade-specific convolutional neural networks

The CNNs from the different network class clades have the same architecture, which is shown in Fig. 2. Each of the CNNs are connected to the input layer. The input layer passes n input examples as a 3D tensor of size n×4000×1 to seven 1D convolutional layers within the different CNNs. The input data are passed as a 3D tensor, as this is required by Keras’ 1D convolutional layers. See the section “Generating data for training and testing” for a more detailed explanation of the input data size.

Open in a separate window

Figure 2:

Architecture of network class-clade convolutional neural networks. All layers with dark-gray backgrounds (except for the input layer) have weights, whereas the layers with white backgrounds are mathematical operations or dropout layers. Below the name of each layer (except for the dropout layers), the sizes of the layers’ output tensors are written. Arrows show the flow of information and how the layers are connected. What constitutes the convolutional neural network (CNN) is encapsulated by a light-gray box. “Conv4” means a 1D convolutional layer with a kernel size of four one-hot-encoded amino acids. The same logic applies to the remaining convolutional layers. “Dense500” means a dense layer with 500 nodes.

The convolutional layers have 40 kernels with no bias parameters and they convolve over 4, 6, 10, 16, 22, 30, and 40 one-hot encoded amino acids, respectively. Their strides are 20, which is the length of 1 one-hot encoded amino acids. The convolutional layers perform valid convolutions. The 3D tensors that are outputted by the convolutional layers are max pooled along their second axes. The max pooled convolutional output values are concatenated, such that a (n×280) tensor is created. Each of the n input examples are then normalized by a division of their summed absolute values. This is done to ensure that the network focuses on amino acid patterns and does not predict based on sequence length. Dropout with a rate = 0.2 is applied on the normalization layer. The normalization layer is connected to a chain of three dense layers with 500 nodes where dropout is applied on each of them (rate = 0.5). See the section “Training settings and network initialization” to read more about the used regularization techniques. In the top of Fig. 2, the CNN is connected to an output layer. This could, for instance, be “CNN1” that is connected to “Output 4_1” in Fig. 1B. All layers with weights that are not output layers consist of Parametric Rectified Linear Units (PReLU) [36].

Generating data for training and testing

Models based on 10-fold cross-validation (CV) were trained. To do this, the peptide classes were divided into 10 bins of similar sizes. To distribute the peptides classes uniformly across the bins, all peptides belonging to all specific sets of classes were counted. If the count value was greater than 10, the peptides were split into the 10 bins, and if the count value was lower than 10, the peptides were uniformly distributed at random across the bins. Like in Reference [11], all peptide sequences were one-hot encoded before introduced to the MultiPep. All sequences were zero-padded until they reached a length of 200 one-hot encoded amino acids.

Training sets, validation sets, and testing sets with a size distribution of 0.8, 0.1, and 0.1, respectively, were generated. For all peptide sequences, labels or targets were created. One indicated membership of a class and zero indicated that they were not members of the specific class.

Training settings and network initialization

The start parameters of the convolutional layers were initialized using orthogonal weights. The initial weights of the dense non-output layers were sampled from a uniform distribution with min=-1 and max=1. The initial weights of the output layers were sampled from a uniform distribution with min=0.001 and max=0.05. The bias parameters of the dense layers and output layers were initialized as zeros.

The Adam [37] update function was used with a learning rate of 0.0005. The remaining hyper-parameters for the Adam update function were set to the Keras default values. The batch size was 128 and dropout with a rate of 20% and 50% on the normalization layer and the dense layers, respectively, was applied (Fig. 2). As a final regularization technique, early stopping after 20 epochs was used. All data were shuffled before an epoch was started. To acquire optimal models, weights for the network class clades were saved whenever their performances on the validation set had improved. Additionally, the weights of the entire network were saved whenever the overall performance had increased. The early-stopping count-down was set to zero whenever the parameters of a “network class clade” were saved. The performance was measured in terms of loss.

The customized MCC loss function

To increase the loss penalty during training, the BCE loss function was extended with a customized version of the MCC. The MCC function is given by

where TP=true positives, TN=true negatives, FP=false positives, and FN=false negatives. MCC values range from −1 to 1, where −1 indicates complete disagreement and 1 indicates a perfect agreement [38]. Here, the MCCloss function for a single sigmoid output node is calculated by

where TP*, TN*,FP*, and FN* are modified under certain circumstances. If yi=0 for all instances in a mini-batch, then TP* is set to be minimum 1. If yi=1 for all instances in a mini-batch, then TN* is set to be minimum 1. If y^i=0 for all instances in a mini-batch, then FP* is set to be minimum 1. If y^i=1 for all instances in a mini-batch, then FN* is set to be minimum one. These modifications are done to avoid that the denominator of the MCCloss function becomes zero.

Again, with a single sigmoid output node, the loss for a single predicted mini-batch ranges between 0 and 2. Zero is given to a perfectly predicted mini-batch and two is given to perfectly miss-predicted mini-batch that contains both yi=1 and yi=0 instances. In the situation where yi=1 and y^i=0 or yi=0 and y^i=1 for all instances in a mini-batch, the loss can never reach exactly 2 due to the TP* and TN* definitions. So, to avoid “losing” loss while training, the loss is set to 2 if |∑iNyi- ∑iNy^i|=N, where N is the number of examples in the mini-batch. In the situation where yi=1 and y^i=1 or yi=0 and y^i=0 for all instances in a mini-batch, the loss can never reach exactly 0 due to the FP* and FN* definitions. To circumvent this, the MCCloss was implemented such that its outputs are multiplied by 0 if |∑iNyi- ∑iNy^i|=0.

The complete loss function with MCCloss

The entire loss function for a single output node given all instances of mini-batch is defined as

where N is the number of examples in the mini-batch. The loss function can be seen as two-step function where the MCCloss values are calculated for the entire batch in the first step and then added to the BCE loss values, which then are averaged. During the updating steps while training, the loss is summed across all output nodes.

Loss function with class weights

Models using class weights were trained, which is a standard way of training models on data with imbalanced class-size distributions [39]. The class weights were based on the data from the training sets and were generated using the “compute_class_weight” function from Scikit-learn [40]. This means that a class weight, Clj, was generated for every j class. While training, the loss of a single output node was calculated using

where N is the number of examples in the mini-batch and Clj is the class weight of the class in question. During the updating steps while training, the loss is summed across all output nodes.

Webtool and stand-alone-program

A Numpy [41] version of MultiPep based on the trained weights was implemented. It is available on https://agbg.shinyapps.io/MultiPep/. Users can choose between using the average of all CV models’ prediction as the final output or to get the max of all CV models’ predictions of the different classes as the final output. This is true for the webtool and the stand-alone-program. The webtool rounds the predictions to three decimals, while the stand-alone-program provides a higher number of decimals. The stand-alone-program is available at https://github.com/scheelelab/MultiPep. If predicting more than 500 in one batch of peptide sequences, we recommend using the stand-alone-program.

MCC loss function rationale

It has been suggested that the MCC function is a good performance metric for imbalanced data sets [42]. Thus, inspired by the work in Abhishek and Hamarneh [43], we found it intriguing to use a customized version of MCC function as parameter-free loss function for our data set with an imbalanced class–size distribution. However, the MCC function as is was in our opinion not suited for a classifier that, while training, iteratively looks at mini-batches of training examples. Especially not when training on data sets with small classes. For example, the MCC function will include division by 0, if for example, only true negatives or true positives are predicted in a mini-batch (Equation 6). Our loss function is a combination of the standard BCE and our MCCloss function, as we found that this combination yielded the best results. In this context, the MCCloss function produces extra penalties for the predictions with a low MCC score.

Model performances

Using data derived from multiple peptide databases, we conducted 10-fold CVs and tested the models on their associated test sets. We trained several models using two different loss functions. The neural network architecture and all hyper-parameters and regularization settings were always the same. First, we trained models using our loss function with the MCCloss penalty, where we applied our training scheme that saved the parameters of the best “network class clades” whenever performance on the validation sets increased. Secondly, while training these models, we saved the parameters of the entire network whenever the overall performance had increased. Thirdly, we ran a 10-fold CV using the weighted BCE loss function. We used the training scheme, where the parameters of the “network class clades” were saved whenever their performances increased and the training scheme where we saved all network parameters whenever the overall performance had increased. This produced 4 × 10 different models, which we used for a comparative study to find the best way of training the network.

Tables 4 and ​and55 show the mean and standard deviations of the performances of the trained models on the CV test sets using MCC and the F1 score as performance metrics. In the tables, the mean values have been compared and the standard deviations are stated to provide readers with an overview of the robustness of the models. Similar tables showing the performance on the test sets, but using precision, recall, and accuracy as performance metrics are available in the Supplementary Information (Supplementary Tables S2–S4). Additionally, plots of receiver operating characteristic (ROC) curves of MultiPep models’ predictions of the test sets can be found in the Supplementary Information (Supplementary Figs S5–S14). Tables showing the performance on the validation sets can as well be found in the Supplementary Information (Supplementary Tables S5–S9). The classification threshold was 0.5 and all mean values in the tables have been rounded to three decimals, whereas all standard deviation values have been rounded to two decimals. We did not use the upper layers’ predictions for the comparisons, as we were only interested in finding how the models performed on the bioactivity classes. The performances on the test sets and validation sets of the individual CV models trained using the “save the individual network class-clade” training scheme can be found in Supplementary Tables S10–S29.

The performance comparisons suggest that the approach where we save the parameters of the network class clades whenever their performances increase is slightly better than saving all parameters when whenever the overall performance has increased. However, if taking the performances on the validation sets into account, it seems to be favorable to use the “save the individual network class clade” approach. Further, the study indicates that using our combined BCE and MCCloss loss function in general produces better classifiers than when using weighted BCE as loss function. Altogether, the performance of our models trained using the “save the individual network class clade” training scheme seems to produce reasonable and robust classifiers with an average MCC score above 0.5, except for the “antiparasite,” “antifungal,” and “antioxidative” class.

Comparisons against state-of-the-art binary classifiers

We next aimed to compare MultiPep against existing state-of-the-art binary bioactivity classifiers. We could not find a range of tools necessary for comparing all bioactivity classes predicted by MultiPep. Thus, we chose a few different classifiers based on their performances and bioactivity focus (Tables 1 and ​and6).6). We only used the tools’ pre-trained models available on the tools’ webpages or on GitHub (Table 6). As many of the chosen classifiers did not take peptides of length 2–200, it was necessary to create different sets of peptide sequences that were compatible with the tools’ input formats (see Supplementary Table S30 for sizes of peptide sets).

For the tool by Veltri et al. [8], Deep-AmPEP30 and RF-AmPEP30, the positive peptides were from the “antibacterial” class whereas the negative peptides were drawn from the classes “cell–cell signaling,” “neuropeptide,” and “peptide hormone.” For mACPred, ToxinPred, HLPpred-Fuse, and AnOxPePred, the positive peptides were from the classes “anticancer,” “toxic,” “hemolytic,” and “antioxidative,” respectively, and the negative peptides were as well drawn from the classes “cell–cell signaling,” “neuropeptide,” and “peptide hormone.” For PredNeuroP, we sampled positive peptides from the “neuropeptide” class and negative peptides from the classes “antiparasite,” “anticancer,” “antibacterial,” “antifungal,” “antimicrobial,” and “antivirus.” For the tool by Veltri et al. [8], we used their newest model for the predictions. For ToxinPred, we used the tool’s Swiss-Prot-based SVM model. For AnOxPePred, we used the tools “Peptide Mode” and only used its “Free Radical Scavenger”-prediction mode.

Additionally, we compared MultiPep with THPep, a tool specialized in finding tumor homing peptides [44]. Here, the positive peptides were from the class “anticancer” and the negative peptides were from the classes “cell–cell signaling,” “neuropeptide,” and “peptide hormone.”

For all peptide sets created, we ensured that there were no overlaps between the positive and negative peptides. We did not make any effort to create balanced data sets, as this should not matter for the comparisons. For the comparisons, we used the MultiPep model trained using our “save the individual network class-clade” approach that produced the lowest overall loss on the test set (Supplementary Tables S15 and S25). All peptides from the above-described peptide sequence sets were drawn from the model’s test set, as the peptides in this set contain unseen peptides that have not been part of model’s training. We did not consider if the sequences in the smaller sequence sets were part of the different tools’ training sets. All peptide sequence sets are available on https://github.com/scheelelab/MultiPep. We observed that MultiPep in general outperforms existing peptide prediction tools when subject to our model’s test set (Table 6 and Supplementary Table S31). While the recall ability of MultiPep is bested by the other tools except for Deep-AmPEP30, RF-AmPEP30, and ToxinPred, MultiPep consistently outperforms the other tools on the MCC score, precision, and accuracy.

Comparing prediction error of MultiPep against PeptideRanker

As an additional analysis, we compared MultiPep with PeptideRanker. PeptideRanker is a binary classifier that can predict general bioactivity of peptides. Thus, we found it interesting to compare MultiPep against this tool. Using PeptideRanker, we predicted all peptides of the test set of the model with the lowest loss on the test trained using our “save the individual network class-clade” training scheme (Supplementary Tables S15 and S25). The set contained 4585 peptides and we did not consider whether any of the used peptides from out test set were part of the training data of PeptideRanker. For both MultiPep and PeptideRanker, instead of using peptides with no bioactivity, we used all peptides of the test set as true positives and found the error of the tool’s predictions. We calculated the error of PeptideRanker’s predictions by ∑i=1N1-roundpN, where N is the number of examples in the test set, p is the predictions, and round() is a function that rounds values to their nearest integer. The rounding function implies that the threshold for the predictions was 0.5. For MultiPep, we took the predictions on the Output_1 output layer and found how they diverged from the labels of that layer. We calculated ∑i=1N∑j=1k|yj-round(pj)|N, where N is the number of examples in the test set, is an element of the row vector y with labels, is an element of the row vector p with predictions, and round() is a function that rounds values to their nearest integer. Also, to make a more direct estimate of whether MultiPep can assign a prediction above 0.5 to bioactive peptides, we calculated an additional error percentage by ∑i=1N1-roundmax(p)N, where N and p are the same as above and max() is a function that finds the maximum value of the row vector p. With this approach, we calculated the prediction errors (Table 7). This analysis shows that MultiPep, with a low prediction error, successfully identifies a higher number of bioactive peptide from a list of peptides with known bioactivity.

MultiPep versus peptipedia

Peptipedia is a tool that can predict many different bioactivity classes simultaneously. The authors describe how they have trained 44 different RF models that each can predict belonging to a specific bioactivity class [27]. As described for the comparisons above, we used the MultiPep model trained using our “save the individual network class-clade” approach that produced the lowest overall loss on the test set (Supplementary Tables S15 and S25). We used the test set of this model to generate a test-subset that was used as input to the Peptipedia models (Supplementary Table S30). Noteworthy, as the Peptipedia webtool was unavailable during the writing of this article, we retrieved the Peptipedia predictions via personal communication with the authors of the article presenting Peptipedia [27]. In the prediction data received by the Peptipedia authors, 41 bioactivity classes were included. The Peptipedia prediction is available on https://github.com/scheelelab/MultiPep. Both MultiPep and Peptipedia can predict many different bioactivity classes. Supplementary Table S32 presents an overview of classes that were deemed comparable and incomparable, which was based on the name of the classes. For this grouping of the classes, it was taken into account how the MultiPep classes were generated (Supplementary Information). If a MultiPep bioactivity class was deemed comparable with more than one Peptipedia bioactivity class, the predictions of these Peptipedia classes were merged such that, for every predicted peptide, the max value of the predictions was used. Comparable bioactivity classes are expected to contain similar and identical peptides. Table 8 shows how MultiPep and Peptipedia performed on the generated test-subset. With the exception of the recall score for the antiparasite and the antifungal classes, we observed that MultiPep outperforms Peptipedia on all performances for predicting the compared bioactivity classes.

MultiPep models versus logistic regression models

We next aimed to compare our MultiPep deep learning models with linear logistic regression models. Therefore, we trained 10 logistic regression models using our peptide data and our 10-fold CV scheme. We used the logistic regression functionalities provided by the Python package Scikit-learn [40]. We initialized the following hyper-parameters as follows: random_state = 1234, max_iter = 1000, and class_weight = “balanced.” The remaining hyper-parameters were set to their default values. The peptides were one-hot encoded before introduced to the logistic regression models. When comparing the performances of the two model types, it can be seen that MultiPep outperforms the logistic regression models (Tables 4 and ​and55 and Supplementary Tables S2–S9, S33, and S34). Nevertheless, the logistic regression models seem to have a generally better recall performance.