iAVPs-ResBi: Identifying antiviral peptides by using deep residual network and bidirectional gated recurrent unit

1.   Introduction

Viruses are microscopic infectious complexes that replicate in host cells. As a type of cellular organism, viruses have a very simple structure and do not have their own metabolic system, which means they have to colonize living cells ^[1]. As a result, once the virus has left the host cell, it loses all its vital activities and can no longer replicate, but it is highly pathogenic ^[2]. The great diversity of viruses, their genetic variation, their unclear routes of transmission and their tendency to acquire specific resistance to antiviral drugs through the evolution of mutations make it extremely difficult for researchers to develop effective and safe specific treatments ^[3,4]. Infectious diseases caused by viral pathogens have a major impact worldwide. Zoonotic viruses such as Ebola, Zika fever, West Nile virus, HIV, SARS-CoV and SARS-CoV-2 are particularly dangerous. COVID-19, in particular, has caused millions of deaths worldwide. The development of antiviral drugs is therefore essential ^[5].

Peptides have interesting pharmacological properties such as high selectivity and relative safety ^[6]. Compared to conventional non-peptide drugs, these potential antiviral agents have a number of advantages. Peptides are highly specific, readily present in the human body and perform a wide range of biological functions. Furthermore, peptides are inexpensive to produce and easy to modify and synthesize ^[7]. They are mainly used as signaling and regulatory molecules in various physiological processes ^[8]. Previously, antiviral peptides (AVPs) are isolated from animal excretions after activation of the host defense process. Today, however, they are now available by reasonable design, either chemically ^[9] or from recombinant libraries. ^[10]. AVPs can be divided into two categories according to their mode of action on the host organism: virus-targeted and peptide-targeted ^[11]. To inactivate specific viral proteins ^[12], virus-targeted peptides inhibit viral enzymes involved in transcription and replication ^[13,14]. Only a few compounds with antiviral activity are known. AVPs are a kind of small polypeptide molecules with biological activity, which can kill or inhibit viruses. And they have been shown experimentally to prevent viruses from attaching to and invading host cells ^[15,16]. The unique molecular structure and mechanism of action have made antiviral peptides a hot spot in antiviral research. The majority of antiviral drugs act on specific regions or components of the virus to inhibit its growth. New antiviral drugs are being discovered by targeting various stages of the virus's life cycle, such as the process by which the virus enters the host and the process by which it is synthesized within the host ^[17].

The limited availability of therapeutic molecules targeting many viral infections means that new antiviral drug candidates need to be found to control resurgent and drug-resistant pathogenic viruses ^[18]. Experimentally validated antiviral peptides can therefore be used as an alternative strategy against medically important viruses ^[19]. Machine learning algorithms, in particular deep learning algorithms, have been used to efficiently identify antiviral peptides. Thakur et al. ^[20] developed the first anti-viral peptide prediction tool, called AVPpred, using amino acid composition and physicochemical properties and support vector machine. Chang et al. ^[21] showed that a physicochemical model using random forests was better at identifying antiviral peptides. Zare et al. ^[22] classified antiviral peptides using pseudo amino acid composition (PseAAC) and Adaboost. Lissabet et al. ^[23] developed AntiVPP 1.0 by employing random forest algorithm to predict antiviral peptides based on net charge, number of hydrogen bond donors, molecular weight and hydrophilicity index. Schaduangrat et al. ^[24] used different machine learning algorithms and proposed Meta-iAVP, which can extract efficient feature representations based on prediction parameters obtained from feature types. Chowdhury et al. ^[25] proposed Firm-AVP based on the physicochemical and structural properties of amino acid sequences. Li et al. ^[26] proposed DeepAVP, a two-channel deep neural network integration model for the analysis of variable lengths of antiviral peptides. Pang et al. proposed a two-stage classification scheme and established PreAntiCoV ^[27] and AVPIden ^[28] models. Timmons et al. ^[29] proposed ENNAVIA, a neural network-based antiviral peptide identification model. Agarwal et al. ^[30] elaborated on the identification and validation of antiviral peptides. Charoenkwan et al. ^[31] summarized machine learning methods for virus detection directly from sequence information. Manavalan et al. ^[32] evaluated highly specific virus predictors. Kurata et al. proposed a method called iACVP ^[33] and used conventional features, binary profiling and word2vec (W2V) embedding to detect coronavirus antigenic peptides (ACVP). Recent studies have collected AVPs from the literature and published them in databases such as AVPdb ^[3], DBAASP ^[34], CAMP ^[35] and APD3 ^[36]. These databases form the basis of computational predictors that help researchers identify AVPs, save time and reduce labor costs. A number of prediction tools have been developed to identify AVPs, but they still have certain limitations: they only support the calculation of peptides in a certain range, only some models use deep learning and the accuracy of the prediction needs to be improved.

In this paper, the iAVP-ResBi model is proposed to identify AVPs. Firstly, to reflect the information of AVPs more comprehensively, we use the k-spaced amino acid pairs (KSAAP), enhanced grouped amino acid composition (EGAAC) based on the N5C5 sequence, encoding based on grouped weight (EBGW), and composition, transition and distribution (CTD) to extract features. Then, the four features groups are fused by bidirectional long short term memory (BiLSTM). After that, the fused features are input into deep algorithm framework to identify AVPs. The deep learning framework is constructed by combining improved residual neural network and bidirectional gated recurrent neural network (BiGRU). Residual neural network is a kind of shortcut connection network that prevents gradients from exploding and disappearing. We adopt two residual blocks, each consisting of batch normalization (BN) layer, and Relu activation function and two convolution layers. After processing by the residual neural network, the feature information is entered into BiGRU. Unnecessary information is deleted and useful information is retained by means of update and reset gates. Then, we add the softmax activation layer for classification. 5-fold cross-validation and independent validation are used for verifying the generalizability of the model. The results indicate that, compared with previous studies, our model achieves optimal results and significantly outperforms existing models, and the iAVP-ResBi model can accurately identify AVPs. In order to better explain the model constructed, we draw a flowchart of the model building process, as shown in Figure 1.

2.   Materials and methods

2.1.   Datasets

In this paper, to facilitate comparison with existing models, we select the well-designed AVP datasets created by Thakur et al. ^[20]. They extract 1245 peptide sequences from more than 80 relevant papers and patents containing more than 30 types of antiviral activity sequences like HIV, HCV, HSV, RSV, SARS-CoV, influenza, etc., which contain duplicate peptides or identical peptides. Nearly 91% of peptides come from nature, and the rest can be obtained by artificial synthesis. After removing the duplicate peptides from 1245 peptide sequences by sequence alignment, 604 high-efficiency AVPs and 452 least or non-effective AVPs are finally obtained, which are randomly divided into training set ${T^{544p + 407n}}$ (544 positive and 407 negative for cross-validation) and test set ${V^{60p + 45n}}$ (60 positive and 45 negative for the independent validation set), respectively, according to the proportion of 90 and 10%.

In an earlier antimicrobial peptide prediction model ^[37], a new dataset is constructed using non-experimental negative peptides instead of experimentally validated negative peptides in order to design a better model and better test the generalizability of the model. Therefore, Thakur et al. ^[20] construct the other two datasets based on non-experimental negative peptides: the training set ${T^{544p + 544n}}$ (544 positive and 544 negative) and independent validation set ${V^{60p + 60n}}$ (60 positive and 60 negative).

In order to further test the generalizability and ensure the reliability of the model, dataset ENNAVIA-C constructed by Timmons and Hewage ^[29] is adopted for identification of anti-coronavirus peptides. ENNAVIA-C includes 109 peptide sequences with anti-coronavirus activity and 356 peptide sequences with experimentally validated poor or no antiviral activity.

2.2.   Feature extraction

Currently, the following feature extraction methods have been used in the existing studies targeting identification of antiviral peptides: amino acid composition (AAC) ^{[20,21,22,23,24]}, pseudo amino acid composition (PseAAC) ^{[25,26,27,28]}, physicochemical property-based feature ^{[27,28,29,30]}, amphiphilic pseudo amino acid composition (Am-PseAAC) ^[24,31], dipeptide composition (DPC) ^[24,25], g-gap dipeptide composition (GDC) ^[24,29], one-hot coding ^[26], composition of k-spaced amino acid group pairs (CKSAAGP) ^[27,28], binary profile (BP) ^[31,32,33], word2vec encoding ^[33] and BLOSUM62 encoding ^[33].

In order to broaden the sources of features information and overcome the uniqueness and one-sidedness of extracted features, it is necessary to delve as deeply as possible into the internal relationships expressed in the data. In general, deep and efficient features can provide a strong basis for effectively improving classification accuracy ^[38]. In this paper, the KSAAP, EGAAC based on the N5C5 strategy, EBGW and CTD are used for extracting features. Compared with K-mer, KSAAP can describe localized fragment sequence information for amino acid pairs separated by k residues. Compared with GAAC, EGAAC can reflect both components and physicochemical properties, while reflecting the localized nature of the protein sequence. It is now well established that the most important influence on protein folding is the uniqueness of amino acid residues, and our choice of EBGW can reflect the distribution of residues with the same unique characteristics, as well as depict the nature of the protein sequence. CTD can account for differences in properties in the presence of different patterns of amino acid distribution. These feature extraction methods have been widely applied to the identification of other therapeutic peptides ^[39,40] and post-translational modification sites ^[41,42,43].

2.2.1.   K-spaced amino acid pairs

K-spaced amino acid pairing (KSAAP) is an efficient feature extraction strategy that can be used to highlight and identify motifs in protein fragments and sequences ^[44,45,46]. KSAAP can be used to identify flexible and rigid regions of proteins and has been successfully applied to identify a variety of post-translational modification sites ^[42,43]. This coding provides a valuable descriptor for the accurate classification of AVPs and non-AVPs on the basis of short-term interactions of residues in the sequence. The detailed procedure of KSAAP is described as follows ^[43,46]. For a peptide sequence fragment, it calculates the frequency of the occurrence of amino acid pairs separated by k residues ^[47]. KSAAP is described as follows:

Here, we choose k = 0, 1 and 2. For example, ${M_{AxxA}}, {M_{AxxC}}, {M_{AxxD}}, \cdots, {M_{YxxY}}$ represent the number for 2-spaced amino acid pairs. ${N_k} = L - k - 1, \; k = 0, 1, 2$, where L represents the length of the peptide sequence. For each value of k, ${F_k}$ represents the features for k-spaced amino acid pairs. All protein sequences studied in this paper are composed of 20 natural amino acids, that is, no virtual amino acids, so the dimension of each feature space is 20 × 20 = 400. ${F_0}$, ${F_1}$ and ${F_2}$ are calculated as shown in Eqs (1)–(3), respectively. ${F_0}$, ${F_1}$ and ${F_2}$ are combined to get the final $400 \times 3{\text{ = }}1200{\text{ - }}$ dimensional feature space F, as shown in Eq (4).

2.2.2.   Enhanced grouped amino acid composition

The enhanced grouped amino acid composition (EGAAC) algorithm is an improved approach to the GAAC algorithm put forward by Chen et al. ^[48]. The EGAAC algorithm converts character information of peptide sequences into digital vectors. It has been successfully used to predict lysine in succinate sites ^[49]. As the peptide sequences used in this study vary in length, with the shortest sequence consisting of six amino acid residues, EGAAC is used to extract isometric features from the N5C5 sequence, because the N and C terminal residues are important for the structure and function of bioactive peptides ^[50,51]. The N5C5 sequence is constructed by truncating five amino acids from the N-terminal and C-terminal parts of the peptide sequence, and combining them into a new sequence with ten amino acid residues. Compared to fill a short sequence with virtual amino acid 'X' to the longest sequence, the N5C5 sequence strategy can avoid feature redundancy. At the same time, this method better preserves useful information, as the N-terminal and C-terminal parts of the sequence contain more essential information.

According to the five physical and chemical properties of amino acids, Lee et al. ^[52] classify 20 types of standard amino acids into the following five categories:

Aliphatic group: ${g_1}$ = {G, A, V, L, M, I}; Aromatic group: ${g_2}$ = {F, Y, W}; Positive charge group: ${g_3}$ = {K, R, H}; Negative charge group: ${g_4}$ = {D, E}; Uncharged group: ${g_5}$ = {S, T, C, P, N, Q}.

EGAAC scans along the N5C5 sequence with a fixed size window n. The calculation formula is as follows:

where $N(g, n)$ represents the number of amino acids belonging to the g-th group in the window n, and $N(n)$ represents the length of the window. The value of window length is from 1 to L, since EGAAC is analyzed on N5C5 sequence, that is L = 10. In this paper, n is set to 5 by default, and we finally obtain a $(L - n + 1) \times 5 = 30$-dimensional vector for each peptide sequence.

2.2.3.   Encoding based on grouped weight

EBGW is a coding scheme ^[53,54] proposed by Zhang et al. ^[55], which characterizes sequences on the basis of the physical and chemical properties of the amino acids. On the basis of their physical and chemical properties, the 20 amino acids are divided into four groups:

Neutral and non-polarity residue group: ${G_1}$ = {A, F, G, I, L, M, P, V, W}; Neutral and polarity group: ${G_2}$ = {C, N, Q, S, T, Y}; Acidic group: ${G_3}$ = {D, E}; Basic group: ${G_4}$ = {H, K, R}.

The above four sets are combined to get three combinations, ${G_1} + {G_2}vs.{G_3} + {G_4}$, ${G_1} + {G_3}vs.{G_2} + {G_4}$, ${G_1} + {G_4}vs.{G_2} + {G_3}$. Given a peptide sequence $P = {p_1}{p_2} \cdots {p_L}$ with length n, it can be transformed into three binary sequences as follows:

where ${p_i}(i = 1, 2, \cdots, L)$ belongs to any one of 20 amino acids.

Then, each sequence is cut into J subsequence fragments. For instance, for ${H_1}$, the feature of $j$-th subsequence is defined as the following:

where $Sum(j)$ is the number of 1 s in the jth subsequence, $D(j) = {{\text{int}}} (j \times L/J)$ refers to the length of the jth subsequence, ${\text{int}}(\cdot)$ is a function that rounds a numerical value down to the nearest integer. The EBGW scheme defines a peptide sequence as the $3 \times J$ dimension vector. Here, J is selected as 1, 2, 3, 4 and 5. Therefore, the dimension of the EBGW-based feature vector is $3 \times 15 = 45$.

2.2.4.   Composition, transition and distribution

The CTD scheme is originally proposed by Dubchak et al. ^[56] to predict protein folding classes. This feature extraction method can be used to describe the overall amino acid composition of each peptide sequence ^[57]. All amino acids are classified into three categories: polar, neutral and hydrophobic. Each amino acid is then coded as 1, 2 or 3, depending on the category to which it belongs. Composition (C) describes the overall percentage composition of 20 natural amino acids, which is defined as:

where ${n_s}$ is the number of the sth class in the sequence and L is the length of the peptide sequence.

The transition (T) represents the percentage frequency of one type of amino acid followed by another type of amino acid, which can be described as follows:

where ${n_{xy}}$ and ${n_{yx}}$ are the number of dipeptides encoded as 'xy' and 'yx', respectively.

Distribution (D) is to calculate the positions of the first, 25, 50, 75 and 100% of each type of 20 natural amino acids, and the descriptor ${E_i}$ is defined as:

where ${P_1}, {P_{25}}, {P_{50}}, {P_{75}}$ and ${P_{100}}$ respectively measure the position of the first residue, and the occurrence rate of X at 25, 50, 75 and 100%, respectively.

According to the seven physical and chemical properties of hydrophobicity, van der Waals volume, polarity, polarizability, charge, secondary structure and solvent proximity ^[58], the dimension of the feature vector based on CTD is finally $(3 + 3 + 3 \times 5) \times 7 = 147$.

2.3.   Feature fusion

Feature fusion is an effective strategy for processing multiple features. It allows redundant information between different features to be eliminated and discrete information to be merged. Here, BiLSTM is used to combine the four feature types mentioned above.

LSTM is based on the principle of recurrent neural network (RNN) ^[59], where each unit is connected in turn to constitute a directional cycle. After this connection, the internal state of the network has been well constructed. At the same time, LSTM can avoid vanishing gradients, exploding gradients and poor ability to rely on long-range information. The LSTM units dynamically adapt to the length of the studied sequence. Each LSTM unit consists of an input gate, a forgetting gate and an output gate ^[60]. The formula can be expressed as follows:

where ${f_t}, \; {i_t}, \; {C_t}, \; {o_t}, \; {h_t}$ represent the forgetting gate, input gate, cell state, output gate and hidden unit state at time t, respectively. ${\tilde C_t}$ denotes the temporary state (candidate cell information). ${W_f}, \; {W_i}, \; {W_C}, \; {W_o}$ are weight matrices, ${b_f}, \; {b_i}, \; {b_C}, \; {b_o}$ are bias vectors, $[\; ]$ means that two vectors are connected and $\odot$ means pointwise multiplication.

In this paper, BiLSTM is a combination of forward LSTM and backward LSTM. BiLSTM performs reverse processing on the input sequence information and recalculates it according to the long-term and short-term memory algorithms ^[61]. Fixed-length time-step outputs in both directions are obtained and the outputs from both directions are combined into a feature vector.

2.4.   Deep learning algorithm

In order to more accurately identify AVPs, we develop a framework of deep learning based on the residual network ^[62] and the bidirectional gate recurrent unit ^[63], which is an improvement on the framework proposed by Zhang et al. ^[64]. As the number of layers in the depth of the network increases, the identification accuracy reaches a saturated state then produces a rapid degradation. And residual neural networks make these layers to fit the residual mapping, rather than having each stacked layer directly fit the desired underlying mapping. He et al. ^[62] have showed that BN and ReLU are full pre-activation, having smaller classification error than the original residual unit. The BN layer in front of the Relu function can bring regularization effect, reduce over-fitting, and thus get higher accuracy. Afterwards we choose the bidirectional gated recurrent neural network to better capture dependencies with large step distances in the sequence. It controls the flow of information through a door that can be learned.

In this paper, first, the normalization process of the batch normalization layers (BN layer) and Relu activation function are carried out. Then, in the convolution processing, each part of the convolution is composed of two convolution layers. In the residual unit, the input and output information are added by means of conv shortcut, and two residual blocks are connected at the same time. The bidirectional gate recurrent unit includes three BiGRU layers, three dropout layers, a flatten layer, two dense layers, and a dense layer as output layer with softmax activation function. The softmax function first converts logits into a probability distribution and then select the node with the highest probability as our prediction class. The deep learning framework is operated in Python 3.8 and TensorFlow 2.5.0 under PyCharm, and the operating system is 64-bit Windows 10.

A residual network is a kind of hop-connected network that skips intermediate layers and passes previous activation values directly to the next network. In this way, the problems of gradient loss and bursting are effectively mitigated and the depth of the network to be learnt is significantly improved. At the same time, the phenomenon of information missing has been basically solved. The general structure of the network is that information is fed sequentially to each layer for processing. The residual unit can directly transfer the input $x$ to the output as the preliminary result by means of shortcut connection, and the output result is $H(x) = F(x) + x$, where $H(x) - x$ is the residual. This means that the cell inputs are added directly to the cell outputs and then activated ^[65], as shown in Figure 2.

The more layers in the network, the more features can be extracted at different levels. However, a high number of layers in the network can saturate and even reduce the accuracy of the training set. The key part of the residual structure shortcut connection is equivalent of a simple identity mapping. In addition, the whole network can still be trained by end-to-end back propagation. Each residual unit can be expressed by the following formulas:

where $F$ is the residual function, which is a stack of two 1-dimensional (1D) convolutional layers with kernel_size = 3. ${W_l} = \{ {W_{l, k}}|1 \leqslant k \leqslant K\}$ is a set of biases associated with the l-th residual unit. The skip connection $h({x_l})$ represents the transform of 1D convolution with kernel_size = 1 and $f$ is the Relu activation function ^[65]. In this paper, we use the improved residual neural network, that is, the BN layer and Relu activation function are calculated first, and then the weight layer is used to calculate. The feature information can be directly transferred from one unit to another, which maintains much of the integrity of the information.

Our deep learning algorithm integrates residual unit and bidirectional gated recurrent unit. Gated recurrent neural network can solve the problem of vanishing gradients ^[66]. It has been widely used to process data. The traditional GRU structure consists of two parts: reset gate and update gate. The most critical step of GRU is to update the memory stage. At this stage, forgetting and remembering occur simultaneously, and update gates are used to control the extent to which information from the previous state is incorporated into the current state ^[67]. Update the GRU cell at each time step t by applying the following equation:

where ${z_t}$ denotes the update gate, ${r_t}$ denotes the reset gate, ${\hat h_t}$ denotes the candidate status, ${x_t}$ is the input vector, ${h_t}$ is the output vector and $\odot$ means pointwise multiplication.

BiGRU is an improvement on the GRU and helps to improve the accuracy of prediction. It consists of two GRUs: a forward GRU model that takes forward input, and a reverse GRU model that learns backward input. The expression is as follows:

where ${\vec h_t}$ is the forward output, ${\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\leftarrow}$}}{h} _t}$ is the backward output, $[\; ]$ means that two vectors are connected and ${y_t}$ is the output of the model.

2.5.   Evaluation metrics

It is not uncommon for overfitting to occur during training, which means the model can match the training data well, but does not give a good result for data outside the training set. If test data is used to adjust the model parameters at this point, the accuracy of the final evaluation results will be affected. It is now common practice to use cross-validation and independent dataset tests to evaluate the training and generalization ability of constructed model ^[68]. In this paper, we choose 5-fold cross-validation and independent validation to estimate the prediction performance of this model. Meanwhile, some valid evaluation indexes are chosen to verify the feasibility of the model. They are sensitivity (Sn), specificity (Sp), accuracy (ACC), precision (Pre), recall (Rec), the receiver operating characteristic curve (ROC), precision-recall curve (PRC) and Matthew's correlation coefficient (MCC) ^{[69,70,71,72]}. They are calculated as follows:

where TP, TN, FP and FN represent the number of true positive, true negative, false positive and false negative, respectively. The ROC curve shows 1-specificity on the x axis and sensitivity on the y axis, and the area under the curve is called auROC (AUC). The PR curve shows the recall on the x axis and the precision on the y axis, and the area under the curve is called auPRC. The area under the two curves ranges from 0 to 1, which is used to evaluate the effectiveness of the model ^[73]. Obviously, the larger the value of auROC and auPRC, the better the performance of the predictor.

3.   Results and discussion

3.1.   Sequence analysis

To analyze which AVPs residues are preferred at which positions, the frequencies of occurrence of the N-terminal and C-terminal residues are examined using two sample sequence logos. Two Sample Logo is a web-based tool that can compare two groups of multi-sequences, show the amino acid composition at a specific position in the sequence, and get their statistically significant differences. The shortest length of AVPs is 6 in this study. Since AVP activity is concentrated in the N-terminal and C-terminal regions of the sequence, we study 1 to 5 N-terminal and C-terminal amino acids by sequence composition preference analysis. We submit the positive and negative samples for 5N-terminal and 5C-terminal amino acids from ${T^{544p + 407n}}$ and ${V^{60p + 45n}}$ datasets to the sample logo online server (http://www.twosamplelogo.org/) ^[74] to generate the two sample sequence logos. The logos are scaled according to their statistical significance threshold of $p < 0.05$ by Welch's t-test. The upper portion (enriched) is represented by positive AVPs, while lower portion (depleted) is represented by negative AVPs (non-AVPs), as shown in Figure 3.

Figure 3(a), (b) expresses two sample logos for the ${T^{544p + 407n}}$ dataset, Figure 3(c), (d) expresses two sample logos for the ${V^{60p + 45n}}$ dataset. From Figure 3, in the ${T^{544p + 407n}}$ dataset, we can see that positive and negative AVPs are significantly different. Positional analysis of 5 N-terminal residues showed that leucine (L) at positions 1 and 4, tryptophan (W) at positions 2 and 3, threonine (T) at position 1, arginine (R) at positions 4 and aspartate (D) at positions 5 are significantly overrepresented compared with other amino acids, while asparagine (N) at positions 2, 3 and 5, valine (V) at positions 3 and 5, alanine (A) and glutamine (Q) at position 4 are significantly underrepresented. In addition, tyrosine (Y) at position 4 is overrepresented, while Y at position 4 is underrepresented, serine (S) at position 1 is overrepresented, while S at position 2 is underrepresented, isoleucine (I) at positions 2, 4 and 5 are overrepresented while I at position 1 is underrepresented and proline (P) at position 5 is overrepresented, while P at position 2 is underrepresented. These results suggest that positive and negative AVPs are significantly different.

3.2.   Identification performance of our model

In this paper, we propose a model named iAVP-ResBi. Firstly, we select KSAAP to extract short-range interaction information between amino acid pairs considering the need to extract features from multiple perspectives. Due to the important information of the sequence being mostly concentrated in the N-terminal and C-terminal, we extract 30 EGAAC features based on N5C5 sequence. The different feature parameters affect not only the size of the feature extraction, but also the subsequent predictive ability of the model. Choosing appropriate parameters is a key step in feature extraction. Selecting a parameter that is too small will lead to incomplete extracted information, while a parameter that is too large may lead to feature redundancy or even dimension disaster. In order to ensure the integrity of information, we choose the step size from 1 to 5 for EBGW, and combine the feature information of five parts. Considering the physicochemical properties of protein, we convert the peptide sequence into a numerical sequence, and extract the features by CTD. Then, the 1422 features are fused by BiLSTM, and we choose the parameter as 400 to get 800-dimensional features. Finally, we design a deep learning algorithm composed of improved residual block and a bidirectional gated recurrent unit for classification, and better experimental results are obtained.

From our model we can see that, the ACC values of the four datasets reaches 95.07, 98.07, 94.29 and 97.50%, which is a significant improvement. At the same time, the values of Sn, Sp, MCC and AUC of the four datasets have been improved to some extent. More specific information of identification results is listed in Table 1. The ROC curve and PR curve for the ${T^{544p + 407n}}$, ${T^{544p + 544n}}$, ${V^{60p + 45n}}$ and ${V^{60p + 60n}}$ datasets are plotted in Figure 4.

3.3.   Analysis of different feature extraction approach

In model recognition, feature extraction is the foundation, and the choice of its method has a vital impact on the prediction effect. In this work, we select four feature extraction methods including KSAAP, EGAAC, EBGW and CTD. These methods extract important information of peptide sequence from different perspectives. The single feature method is compared with the combined feature method and the result is shown in Table 2. Clearly, the combined features show better performance than any other single feature. Thus, the information obtained from different methods helps to improve the predictive power of the model.

3.4.   Analysis of feature fusion approach

The main purpose of feature fusion is to combine features from multiple sources into a better feature representation to improve model performance. We choose to use the deep learning algorithm BiLSTM to fuse the original features. Through BiLSTM processing, the features we obtained have been improved in importance and relevance. To illustrate the necessity of feature fusion, we compare the features before and after fusion, the experimental results are shown in Figure 5. Evidently, the ACC of the four datasets is improved after using the feature fusion method. This means that the classification effect of the whole model is greatly improved. Therefore, it is essential to use the feature fusion method to process features.

3.5.   Performance comparison with different classifiers

The most important part of the binary classification problem is the construction of the classifier, which directly affects the final test result. In this paper, we design a deep learning algorithm named ResNet + BiGRU. The combination of improved residual neural network (ResNet) and BiGRU effectively improves the identification accuracy of the model. Our ResNet used is pre-activation for the full BN layer and Relu activation function based on the original residual network.

The results show that ResNet + BiGRU has better performance than other classification algorithms. To describe the superiority of our deep learning algorithm, we compare ResNet+BiGRU with different classifiers. More details are shown in Tables 3 and 4. In general, the deep learning algorithm has better predictive performance than the machine learning algorithm on the same dataset. In the machine learning approaches, we compare SVM, XGboost, GaussianNB, logistic regression (LR), ANN with ResNet + BiGRU. The ACC values are 80.86, 80.44, 74.24, 76.97 and 80.33% in the ${T^{544p + 407n}}$ dataset, respectively. For ${V^{60p + 45n}}$, the ACC values are 60.95, 61.90, 73.33, 70.48 and 71.43%, respectively. In the deep learning approaches, we compare DNN, CNN, BiGRU, original ResNet, ResNet with ResNet+BiGRU, with ACC values of 84.05, 85.11, 87.21, 94.02 and 94.54% in the ${T^{544p + 407n}}$ dataset, respectively, and with ACC values of 85.71, 88.57, 85.71, 88.57 and 89.52% in ${V^{60p + 45n}}$, respectively. Compared with other deep learning algorithms, ResNet+BiGRU is an obvious improvement. In addition, the experimental results also indicate that the accuracies of the original ResNet are lower than those of ResNet with pre-activation for the full BN layer and Relu activation function for four datasets, respectively.

At the same time, to demonstrate the superiority of our combination of the ResNet and BiGRU, we also choose to combine the ResNet with BiLSTM, and the experimental results are shown in Tables 3 and 4. The final identification accuracy does not reach a more satisfactory result. Therefore, we adopt the ResNet with BiGRU to improve the accuracy. The experimental results show that our classifier has better identification performance and generalization ability. The parameters of different classifiers are shown in Table 5.

3.6.   Performance comparison with different models

We summarize several currently available models for identifying AVPs. To illustrate more visually the advantages of our constructed model iAVP-ResBi, seven models, namely AVPpred ^[20], the model proposed by Chang et al. ^[21], AntiVPP 1.0 ^[23], Meta-iAVP ^[24], Firm-AVP ^[25], DeepAVP ^[26] and ENNAVIA ^[29], are compared using the same training dataset and test dataset. The values of four evaluation metrics for seven models are shown in Table 6.

The values in Table 6 show the predicted results of our method and the existing methods on the four datasets. Clearly, the results of our model have been greatly improved in both training set and test set. The ACC on the four datasets improved by 3.82, 2.17, 0.41 and 1.85% compared with the best model, respectively, which fully demonstrates the power of our method. The results show that our proposed iAVP-ResBi model has better predictive performance than the previous model.

3.7.   Universal validation of the model

In order to further illustrate the generalizability of our constructed model, we use the dataset ENNAVIA-C created by Timmons and Hewage ^[29] to identify anti-coronavirus peptides. ENNAVIA-C consists of 109 peptide sequences with anti-coronavirus activity and 356 peptide sequences that have been validated to have poor or no antiviral activity. As shown in Table 7, the identification accuracy of anti-coronavirus reaches 96.99% in our constructed iAVP-ResBi model, which exceeds several currently available models. The results show that the model we constructed is an effective recognition model with good generalization ability.

4.   Conclusions

As a long-term focus of researchers, antimicrobial peptides have shown great advantages in the fields of medicine and life sciences. Among them, antiviral peptides (AVPs) are an important component of antimicrobial peptides, which can not only cooperate with the components of the immune system, but also have a good therapeutic potential in drug resistance. Currently antiviral peptides are selected from databases such as HIPdb, APD3, CAMPR3 and LAMP, which are experimentally determined and time-consuming. Meanwhile, antiviral peptides have been favored by biopharmacologists in recent years due to their high pharmacological activity, strong targeting, low toxicity and more mature production technology. With the rapid development of artificial intelligence and deep learning algorithms, a large number of peptide vaccines, peptide nutraceuticals, peptide drugs and reagents will be developed and enter clinical trials in the coming years. So it is very necessary to identify AVPs accurately and efficiently. In this study, we develop a novel model called iAVP-ResBi for identification of AVPs. A single method to extract features is not comprehensive. Here we use the method of combining four features: KSAAP, EGAAC, EBGW and CTD. Then to develop a better classifier, after feature fusion through BiLSTM, the residual neural network combined with BiGRU was used to identify AVPs.

The final data shows that the model we constructed has good results. The value of ACC is 95.07, 98.07, 94.29 and 97.50% on the ${T^{544p + 407n}}$, ${T^{544p + 544n}}$, ${V^{60p + 45n}}$ and ${V^{60p + 60n}}$ datasets, respectively. As a result, the iAVP-ResBi model constructed in this paper has certain reference significance. In our future work, we will be committed to establishing the user-friendly and public web-server for our iAVP-ResBi model for the convenience of the wider research community to use. The datasets and codes are freely available at https://github.com/yunyunliang88/iAVPs-ResBi.

Our model is applicable in the identification of bioactive peptides but may not be particularly applicable in the identification of post-translational modification sites of protein or functional prediction of DNA. In addition, if the number of samples exceeds tens of thousands, our model will have some uncertainty, and at the same time, our hardware equipment may not be able to support it, so we need a high-performance computing workstation for processing.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No.12101480), and the Fundamental Research Funds for the Central Universities (No. QTZX23002).

Conflict of interest

The authors declare there is no conflict of interest.