Innovative Alignment-Based Method for Antiviral Peptide Prediction

1. Introduction

Antiviral peptides (AVPs) represent a promising, innovative and unconventional approach in the ongoing fight against viral infections, a persistent global health concern [1]. The ever-present threat of viral outbreak underscores the critical need for effective antiviral treatments. AVPs have emerged as a promising alternative due to their unique ability to target various stages of viral infections [2].

Viruses, notorious for their vast genetic diversity and adept replication within host cells, present immense challenges in disease containment [3]. Traditional antiviral drugs, while effective in certain instances, often have limited reach and may become obsolete due to resistance [4]. This treatment void has catalysed research into AVPs as alternative antiviral agents. AVPs have shown efficacy against numerous viruses, including lethal ones like HIV [5], influenza [6], hepatitis [7], and emerging zoonotic viruses such as Ebola [8] and Zika [9].

The therapeutic allure of AVPs is rooted in their distinctive mechanisms of action. They hinder viral attachment, fusion, and replication, offering a comprehensive strategy against viral infections [10]. Derived from synthetic libraries or sections of natural proteins, these peptides boast characteristics vital for their antimicrobial actions. Their low toxicity, high specificity, and efficiency marking them as potent contenders in medical applications [11].

Discovering and validating AVPs through lab experiments is resource-intensive, but computational strategies are gaining traction for their potential in identifying peptides with antimicrobial properties. Crafting robust computational models is imperative for effectively pinpointing potential AVPs. Currently, the strategy for processing expansive peptide data banks relies on machine learning (ML), enabling in-depth multidimensional data analysis. Conventional ML algorithms, such as SVM [12,13,14], kNN [15], RF [16,17,18], NN [19], and deep learning (DL) algorithms [20,21,22], have shown prowess in discerning patterns in peptide sequences and assessing new ones.

However, a point of contention in this domain is the efficacy of DL models in forecasting AVPs. Many DL techniques necessitate vast experimentally validated peptide sequence datasets, which are often lacking. One remedy is "data augmentation," although it is largely untapped [21]. As García-Jacas et al. highlighted [23], DL does not significantly outpace traditional ML, and their algorithmic outputs frequently intersect. Additional challenges include the overrepresentation of specific sequences and imbalanced data, which can distort evaluations based purely on accuracy. Moreover, reproducibility challenges persist, as not all scientists disclose their source code or datasets, obstructing these methods' widespread acceptance. Therefore, while ML-centric methods hold promise in predicting active peptides, their refinement is an ongoing endeavour [24,25].

Recent studies have showcased a non-trained supervised technique, the Multi-Query Similarity Search (MQSS), for predicting peptide bioactivities, including hemolysis [26], tumor-homing [27], and antiparasitic [28] activities, with impressive results. This method trumps conventional ML methods in several ways: it is user-friendly, does not rely on web server availability, consumes fewer computational resources, and processes sequences with non-standard amino acids or varying lengths. Remarkably, MQSS models (MQSSMs) function without extensive training, relying instead on fine-tuning certain parameters like sequence alignment type and the similarity cutoff value. They can be developed without needing a negative dataset, a significant advantage given the scarcity of validated negative sequences, ensuring the learning phase is not skewed by data imbalances.

In prior research, we explored the AVPs Chemical Space using interactive mining and complex networks via the StarPep toolbox ([29], https://github.com/Grupo-Medicina-Molecular-y-Traslacional/StarPep). We generated, disclosed, and made public diverse reducts of the AVPs Chemical Space by implementing scaffold extraction to Half-Space Proximal Networks formed from the predefined Chemical Space. These scaffolds now serve as a foundation for designing reference/query datasets for MQSSM creation [30]. As the MQSS technique has not been applied to AVPs, our primary objective is to create MQSSMs for AVPs that match or surpass existing predictors, while addressing the shortcomings of ML predictors. The MQSSMs were crafted using the available resources in StarPep toolbox. Our overarching goal is to enhance the methodologies for identifying and advancing AVPs, potentially revolutionizing antiviral therapies. Additionally, we have assembled the most extensive dataset of antiviral sequences to date, essential for developing effective AVP predictors and ensuring robust model validation.

2. Materials and Methods

2.1. The Multi-Query Similarity Search Model (MQSSM). The Overall Approach

The models proposed here are based on similarity searching. Multiple positive sequences are employed as queries to predict the antiviral activity of a target dataset. This one-class model identifies potential positive sequences based on their similarity to the provided references. The construction of these models involves tuning three different parameters:

Query Dataset: This is defined as the reference dataset that the model uses for predictions.

Pairwise Sequence Alignment Algorithm: Smith-Waterman (local alignment) and Needleman-Wunsch (global alignment) algorithms are employed, using the Blossum-62 substitution matrix for calculating pairwise similarity scores.

Similarity Threshold: The MQSS is a fusion model, specifically a group fusion model. In a group fusion, a similarity is computed for each reference sequence q and target sequence t, denoted as S (q, t). All pairwise similarity scores are combined for each sequence t, and the MAX-SIM rule is applied. The similarity scores are then ranked in decreasing order. A Similarity Threshold (c) is set to determine which sequences are considered to have positive activity. For a better understanding of the process, a graphical representation is provided in Scheme 1.

2.2. Construction of Query/Reference Datasets

When constructing a Query Dataset, a critical consideration is its ability to encompass a significant portion of the antiviral active chemical space without disproportionately representing any category of AVPs or neglecting rare sequences. To achieve this, we leveraged the reductions obtained in a previous study [30].

2.2.1. Recalling the Scaffold Extraction Procedure

To ensure a comprehensive coverage of antiviral active chemical space, we established a Half-Space Proximal Network (HSPN) using the 4,663 antiviral sequences stored in StarPepDB [31], recognized as one of the most extensive databases of biologically active peptides to date (36). Several scaffold (representative AVPs subsets) extractions were executed using the resources available in the StarPep toolbox. During this process, we adjusted the alignment algorithm type, the pairwise identity percentage, and the centrality measure (Community Hub-Bridge [32] and Harmonica [33,34]). These adjustments yielded a total of 20 scaffolds or representative AVPs subsets. Subsequently, we thoroughly examined these scaffolds using a Dover Analyzer [35] to assess the extent of similarity overlap between them.

Based on these observations, we selected 10 scaffolds for their representativeness and diversity. Furthermore, an additional five scaffolds were crafted by combining some of the 20 previously generated scaffolds, considering their pairwise identity percentages. These 15 scaffolds served as the foundational datasets used as references for the MQSSMs. The particularities of each of the mentioned scaffolds are briefly explained in Table SI1.1 of the Supporting Information (SI) file 1 (FileSI1).

2.3. Target Datasets for Calibration and Validation of MQSSMs

We conducted a comprehensive literature review to gather a range of positive and negative sequences for testing our models. Ultimately, we selected 15 datasets, containing positive and negative sequences from synthetic and natural origin. These datasets from literature played a crucial role in various stages of model calibration and validation. All the mentioned datasets are depicted in Table 1.

Furthermore, we created two new datasets, amassing a total of 20,136 positive sequences and 70,130 negative sequences. These 90,266 sequences underwent several filters to eliminate redundancy, resulting in an "Expanded" Dataset comprising 55,822 sequences, including 3,178 positive and 52,644 negative sequences. The different sources for the “Expanded” Dataset are specified in Scheme 2.

Additionally, for benchmarking against state-of-the-art methods, we curated a "Reduced" Dataset with specific criteria, excluding non-standard amino acids and restricting sequence lengths to between 7 and 30 amino acids, as commonly specified by many existing predictors. Refer to Scheme 2 for a better depiction of the distinct filters that this dataset went through. The "Reduced" Dataset consists of 27,692 sequences, encompassing 1,419 positive and 26,273 negative sequences. Both, the “Expanded and Reduced” versions can be found at the supplementary materials (File SI2).

2.4. Construction, Selection and, Improvement of MQSSMs

As mentioned earlier, three primary parameters require fine-tuning for the construction of Multi-Query Similarity Search Models (MQSSMs). The first parameter is the choice of the Query Dataset, for which 15 different scaffolds were evaluated. The second parameter is the selection of the alignment algorithm, alternating between global and local alignments. Lastly, the similarity threshold was varied from 0.3 to 0.9. These parameters collectively allowed for the potential creation of 210 variations of MQSSMs. Consequently, reducing this number of models was imperative.

To select the best-performing models, a two-phase approach was implemented, as depicted in Scheme 3. The first phase, known as the calibration phase, aimed to significantly reduce the number of models and identify key trends in parameter calibration. The calibration phase was further divided into two rounds. In the first round, 210 models were assessed using the datasets TS_StarPep, TR_StarPep, and Ex_StarPep, resulting in a reduction to 80 models. The second round of the calibration phase involved testing the models against six different datasets: AVPIden, AMPfun, ENNAVIA-A, ENNAVIA-B, Imb, and Thakur. After this evaluation, a set of 50 models remained.

These 50 models subsequently entered the validation stage, where the primary objective was to fine-tune parameter values while also assessing the models against more specific datasets. Once again, the validation phase was subdivided into two rounds. In the first round, the 50 models were tested against five datasets: Sharma, AI4AVP, ENNAVIA-C, ENNAVIA-D, and Imb_CoV, with the last three containing Anti-SARS-CoV sequences. This round led to the selection of 32 models, which were further evaluated against the Expanded Datasets. Following this analysis, a final set of 12 models was chosen, comprising six that utilized global alignment and six that employed local alignment (M1-M12).

As 12 models still represented a considerable number, an examination was conducted to assess the degree of overlap in the sequences recovered by these different models. After a cursory examination, three base models were selected referred as M3, M7, M12.

2.5. Scaffold Fusion

To enhance the model's performance, a scaffold fusion was performed with the aim of expanding the reference dataset. This fusion process involved combining the individual scaffolds from the best-performing base models (M3, M7, M12) to create a new query dataset. Once this new query dataset was constructed, 14 MQSSMs were generated using it, including variations in the alignment algorithm and the similarity threshold during the creation of these models.

Subsequently, these models underwent the previously described workflow of calibration and validation. Following this process, the best-performing model was identified and is now referred to as M13 (Scheme 3).

2.6. Scaffold Enrichment

To further enhance the base models, a query enrichment strategy was implemented. This involved gathering the positive sequences from the datasets used in Calibration Phase, round 2 (Thakur, ENNAVIA-A, AMPfun, and Imb) into a single dataset. This process resulted in 2403 unique sequences. To ensure the absence of these sequences from the currently top-performing model’s scaffolds, a pairwise similarity comparison was conducted using Dover Analyzer (REF). Subsequently, these sequences were used to construct a Half-Space Proximal Network (HSPN), similarly to what was done for the StarPepDB AVPs sequences.

With the newly curated and validated dataset in hand, it was integrated into the StarPep toolbox to construct an HSPN, similarly to the HSPNs developed previously for the entire AVP space. Following HSPN construction, Scaffold Extraction was applied to the network, varying the centrality measure between Harmonica and Community Hub-Bridge. The alignment algorithm type was adjusted between local and global, considering only 80% and 90% as the sequence identity percentages. This process yielded 8 new scaffolds, designated as” external scaffolds.” These scaffolds were utilized to develop new MQSSMs employing the same methodology as before, resulting in 112 new models. These newly generated models underwent testing against the Calibration and Validation workflow of 15 datasets to identify the most effective external scaffolds.

Finally, two models, named E1 and E2, were selected. The scaffolds from models E1 and E2 were added to those of models M3, M7, M12, and M13, ensuring that there were no repeated sequences in the enriched scaffolds. Once the enriched scaffolds were obtained, the models retained the alignment type and similarity cutoff of the base models, and they were tested against the 14 individual datasets and the expanded dataset (Scheme 3).

2.7. Performance Evaluation

The performance of all the models and state-of-the-art models was measure using the metrics of Accuracy (ACC), Sensitivity (SN), Specificity (SP), Mathews Correlation Coefficient (MCC) , False Positive Rate (FPR) and F1 Score [41]:

$$ACC={\displaystyle \frac{TP+TN}{TP+TN+FP+FN}}$$

$$SN\left(Recall\right)={\displaystyle \frac{TP}{TP+FN}}$$

$$SP={\displaystyle \frac{TN}{TN+FP}}$$

$$\left(MCC\right)={\displaystyle \frac{TP\times TN-FP\times FN}{\sqrt{(TP+FP)\times (TP+FN)\times (TN+FP)\times (TN+FN)}}}$$

$$Precision={\displaystyle \frac{TP}{TP+FN}}$$

$$F1={\displaystyle \frac{2\times Precision\times Recall}{Precision+Recall}}$$

$$FPR={\displaystyle \frac{FP}{FP+TN}}$$

where, TP are the True Positives, TN are the True Negatives, the FP are the False Positives and FN are the False Negatives.

Additionally, the Friedman test was employed to analyse the performance rankings of the various MQSS models across all described metrics. The corresponding significance test was carried out using the KEEL software ([42], https://sci2s.ugr.es/keel/development.phpx1-2000) and the Non-parametric Statistical Analysis Module. This non-parametric test is well-suited for comparing multiple models across different datasets or conditions, especially when dealing with ranked data, as in the present case.

2.8. Comparison with State-of-the-Art

To assess the robustness of the MQSSMs, a selection of state-of-the-art predictors was made, primarily considering their availability through web services and ease of implementation. The Table 2 provides a summary of the predictors used in this section. For the comparison with the predictors the Reduced Dataset was employed as previously stated

3. Results and Discussion

3.1. Performance of MQSSMs at the Calibration Phase

This segment focuses on managing the initial amount of 210 models via a reduction process. In the preliminary phase of calibration, an evaluation of these 210 models was conducted utilizing three datasets: TR Starpep, TS StarPep, and EX Starpep. These datasets, predominantly derived from StarPepDB, demonstrated a relatively uniform behavior across models’ performance, as depicted in Figure 1. Such an outcome is expected considering the considerable overlap in sequences shared between these datasets and the employed "Query."

During the next phase of calibration, six datasets not associated with StarPepDB were incorporated (Figure 1). Predictably, in this phase, the models reduced their effectiveness compared to the initial round. A critical revelation from this stage was a better performance of models when dealing with datasets composed of both randomly generated negative sequences and sequences yet to be experimentally tested. In contrast, when the datasets had experimentally validated sequences, the performance metrics dropped. This behaviour might be linked to the fact that a significant number of these negative sequences closely resembled the positive ones. Consequently, the methods based on alignment faced challenges in differentiating unique attributes of each category.

This problem was especially pronounced in datasets like ENNAVIA-A and Thakur, where experimental sequences were used as negative datasets. Conversely, with ENNAVIA-B, which contained the same positive sequences, the models more effectively identified non-antiviral sequences. It is important to emphasize that experimentally validated negative sequences are scarce in comparison to positive sequences. Hence, the predominant challenge in modelling lay in improving the detection rate (recall) of positive sequences overall.

From the calibration stage, clear trends in the MQSSMs parameters emerged. A notable observation was that models with larger scaffolds performed better. This improvement was attributed to the more detailed characterization of the AVPs' chemical space. Scaffolds such as Md4, Md5, SG4, SG5, SL4, and SL5 generated the most model variants (Table SI1.1). In terms of alignment, global alignment proved more effective with lower sequence identity percentages, while local alignments were more successful with higher percentages. For simpler scaffolds, global alignments consistently outperformed local alignments, regardless of the identity percentage.

3.2. Performance of MQSSMs at the Validation Phase

In the validation stage, we assessed the models against datasets with sequences targeting specific virus like SARS-CoV. Datasets including ENNAVIA-C, ENNAVIA-D, and Imb_CoV were used for this purpose. As shown in the Figure 1, the base models had poor performance with ENNAVIA-C and Imb_CoV. This result was expected because the model references are from 2019. This highlights the importance of the "Query" dataset's representation in model performance. However, the models showed better performance on the ENNAVIA-D dataset, which contained random negative sequences. This aligns with previous observations of model behaviour in similar datasets.

During the model selection process, 32 models chosen in the first round of the validation were tested against the Expanded dataset. From this, 12 models were identified as top performers using a multi-variable Friedman ranking method. These models, labelled M1 to M12, included 6 based on global alignment and 6 on local alignment strategies. The parameters for these models are summarized in the Table SI1.2. Additionally, from these 12 models, the top 3 performing ones were further singled out. This selection focused efforts on fine-tuning these models, particularly towards enhancing the recovery of positive (AVP) sequences.

3.3. Improving MQSSMs Performance by Fusing Scaffolds

Continuing from the focused analysis of the top models for positive sequence recovery, an initial approach to enhancement involved combining the scaffolds from models M3, M7, and M12 (md4, SL5, SG5) into a single, consolidated scaffold. This step included the removal of duplicate sequences, culminating in a scaffold that contained 3206 unique sequences. Subsequent testing of this modified scaffold indicated a slight improvement in the performance of models using global alignment with a 90% similarity threshold. Labelled as M13, a new model was crafted using these parameters. The improvements in M13 were primarily due to the expanded representation of sequence space, enriched by the increased number of sequences.

Although the modifications provided insights into the nuances of the MQSSMs, their overall efficacy remained unsatisfactory. The subsequent strategy aimed to leverage the performance information gathered during the calibration phase. As shown in Figure 1, the base models struggled with datasets such as Thakur, ENNAVIA-A, AMPfun, and AVPiden. This struggle was likely due to the inadequate representation of diverse sequences from these datasets in the MQSS scaffolds. Furthermore, the presence of many experimentally validated negative sequences in some datasets increased the difficulty of making accurate predictions.

3.4. Improving MQSSMs Performance by Enriching the Best Scaffolds

In tackling the identified issues, a Half Space Proximal Network (HSPN) was made by pooling together positive sequences from the challenging datasets (Thakur, ENNAVIA-A, AMPfun, and AVPiden). A total of 2403 sequences were employed in constructing the HSPN. This effort resulted in the production of 8 scaffolds, out of which the top two were chosen to enhance the existing scaffolds. Consequently, this led to the introduction of new, improved models: M3+, M7+, M12+, and M13+. The "+" in their names signifies their references enrichment. Post-enhancement, these scaffolds contained 3155, 3437, 3472, and 3606 sequences, respectively. To incorporate unique scaffolds do not present in StarPepDB, we crafted E1 and E2 models using an external scaffold, comprising 1517 and 1261 sequences, respectively. This increased our analysis pool to 10 models, adding 6 newly enriched models to the pre-existing M3, M7, M12, and M13. Complete details of these models are outlined in Table SI1.2.

Following their development, the 10 models were extensively tested across 15 databases, the 14 datasets from our workflow plus the Expanded Dataset (File SI3). We used a Friedman ranking system to halve the number of top-performing models, evaluating them based on accuracy (ACC), specificity (SP), sensitivity (SN), Matthew's correlation coefficient (MCC), and F1 score. The ranking results identified M3+, M13+, M7, M12, and E1 as the most effective, as highlighted in grey in Table 3.

Notably, while models with a greater number of reference sequences like M3+ and M13+ ranked high, models such as E1, M7, and M12 with fewer sequences performed comparably well. This suggests that the effectiveness of the references hinges more on their diversity and representational range than on their quantity.

3.5. Benchmarking the Best MQSSMs against State-of-the-Art Predictors

The top 5 models were subsequently benchmarked against existing predictors in the literature, providing a comparative assessment of their performance relative to other available tools. To ensure unbiased comparisons, any sequences shared between the scaffolds of models M3+, M7, M12, M13+, E1, and the Reduced dataset were eliminated. This action reduced the number of positive sequences to 116, while the number of negative sequences remained unchanged. It is important to note that many negative sequences in the Reduced dataset are part of the training sets for the external predictors, but these were not removed. This decision placed our models at a significant disadvantage in terms of performance evaluation. We tested a total of 14 external predictors (File SI4), using evaluation metrics such as ACC, SP, SN, MCC, FPR and F1 Score. MCC, which is unaffected by the dataset's imbalance, was the primary metric for our analysis.

The alteration made to the Reduced dataset resulted in a marked decrease in the performance of the MQSSMs. This decline was especially pronounced in the SN and MCC values. However, ACC and SP remained relatively stable, a situation due to the substantial imbalance between positive and negative cases in the class distribution. The notable drop in sensitivity underscores a significant, consistent shortfall in correctly identifying positive sequences. This inability to accurately recall true positives notably affected the MCC, with a stark decrease observed, for example, in the M13+ model where MCC plummeted from 0.731 to 0.214, as detailed in Table 4. The F1 score, which depends on both recall and precision, also experienced a corresponding decline.

Despite the less-than-ideal results, the MQSSMs still surpassed the external predictors in overall performance. Analysing the performance data of the external predictors, Figure 2 highlights two distinct trends. Some models are highly effective at identifying most positive sequences, resulting in a high SN but with the trade-off of a higher rate of false positives. Conversely, other models excel at correctly classifying all negative sequences but tend to misclassify many positive ones, a trait observed in the MQSSMs we developed. Typically, most models based on deep learning fall into the former category, demonstrating high sensitivity, whereas traditional ML models are more likely to be in the latter category, with a stronger emphasis on specificity.

The analysis revealed that no single predictor excelled in all evaluation categories, supporting Garcia-Jacas et al.’s conclusion that DL methods may not be the most effective for AVP prediction [36]. Extending this observation, it's evident that none of the ML models tested demonstrated satisfactory performance, suggesting a need for significant improvements. A central issue identified is the quality and representativeness of the training data. Most positive sequences used in training exhibit high similarity, often up to 90%. Moreover, the experimentally validated negative sequences often mirror their positive counterparts. The challenge, therefore, does not lie in the complexity of the models’ architectures but rather in the data's availability and diversity. This highlights an ongoing challenge in gathering and effectively utilizing comprehensive data for such models.

Evaluating state-of-the-art predictors revealed a common issue of accessibility. Many predictors were difficult to assess for various reasons. A significant issue was the poor construction of several web servers, leading to operational failures or frequent malfunctions. This problem affected even relatively new servers, those less than 2 years old. Furthermore, a few research repositories lacked comprehensive instructions, complicating their use. This difficulty in accessing and implementing these tools echoes concerns previously noted by researchers [46] about the availability of source codes. Nonetheless a comprehensive summary of the prediction tools currently available is provided in the supplementary information (Table SI1.4).

In contrast to these challenges, the MQSSMs distinguish themselves with their accessibility. They are available through the StarPep toolbox standalone software, which boasts a user-friendly interface, making these models more approachable and easier to use for users.

Despite the need for improvements to address their shortcomings, the MQSSMs still have a performance edge over traditional ML models. This is substantiated by a Friedman MCC, ACC, SP, SN, and the F1 Score. An important aspect to note is that the MQSSMs require considerably fewer computational resources and are not constrained by sequence length limitations, presenting significant benefits.

When choosing a prediction model, it’s essential to align with the specific needs of the researcher. Some may prioritize identifying a larger number of potential AVPs, while others might prefer a smaller, more accurate set to reduce false positives. Considering the resource-intensive nature of experimental procedures, the latter approach is often more practical for synthesizing potential AVPs, as it optimizes the balance between resource use and the likelihood of accurate predictions.

4. Conclusions

This research marks a significant advancement with the development of Multi-Query Similarity Search Models (MQSSMs). These models, devised from a deep understanding gained in the chemical space exploration and scaffold extraction phases, employ a novel approach that leverages structural similarities for predicting biological activities. This has been instrumental in improving model selection.

Another major achievement of this study is the assembly of the most extensive datasets of AVPs, referred to as the “Expanded” and “Reduced” datasets. These datasets to were designed to cater to a spectrum of research needs in the field AVP modelling and prediction and are available at the File SI2 in supplementary materials.

Such rich data diversity is essential for the development of effective AVP predictors. Following comprehensive evaluations and filtering processes, five MQSSMs were selected for their superior performance. These models were rigorously tested across various metrics, with the top model demonstrating impressive results: ACC=0.969, SP=0.979, SN=0.802, MCC=0.71, FPR=0.021, F1=0.746. In contrast, 14 contemporary ML-based predictors, though extensively promoted, were surpassed by the MQSSMs. This not only emphasizes the limitations of existing methods but also highlights the advanced capabilities of MQSSMs in predicting AVP sequences, effectively handling the challenges of variable-length sequences and imbalanced datasets. As these models notably eclipsed existing state-of-the-art predictors, they have set a new standard in the field of AVPs discovery.