Peptide Hemolytic Activity Analysis using Visual Data Mining of Similarity-based Complex Networks

1. Introduction

Peptides are relatively small chains of amino acids (AAs) that can be chemically synthesized or purified from living organisms [1]. Our own bodies naturally produce peptides that carry out several critical physiological functions including healing, defense against infections or as chemical messengers [2,3]. Currently, peptides are becoming highly relevant in medical applications as they have shown to exhibit not only promising therapeutic activities such as antimicrobial, antifungal, antiviral, antiparasitic and anticancer but also due to their interesting pharmacological characteristics such as high efficacy, target selectivity and good tolerability [4,5,6]. Peptide drugs were reported to have sales of more than $70 billion in 2019 [7] and in the last decades, they have gained more attention as potential therapeutic drugs than antibodies and small-molecule-based drugs [8,9].

Diseases such as fibrosis, asthma and cancer are treated using peptide-based therapies [2,3]. For instance, the synthetic peptide Leuprolide has been successfully used to treat prostate and breast cancers by acting as an agonist of the gonadotropin-releasing hormone (GnRH) [10]. In addition to Leuprolide, just over 60 therapeutic peptides are commercially available and about 150 are being tested in clinical trials [1,6,11]; however, these numbers are quite low compared with the several thousand potential therapeutic peptides that have been identified [12]. This concerning low proportion of peptide drugs on the market is partially explained by the short half-life, lability during storage, poor oral bioavailability and undesirable toxicity that peptides usually have [2,6,9]. Mainly, peptide-associated hemolysis is perhaps one of the main drawbacks of these potential therapeutic drugs [4] since the products released after the lysis of red blood cells (RBCs) can lead to systemic inflammation and widespread tissue damage [13].

Currently, there are many datasets available containing information about hemolytic peptides. The main databases include: (i) Hemolytik [14], with more than 2000 experimentally validated hemolytic peptides; (ii) Database of Antimicrobial Activity and Structure of Peptides (DBAASP v3) [15], with more than 11321 entries showing information on hemolytic and cytotoxic activities of antimicrobial peptides (AMPs); and (iii) StarPepDB [16] which is a graph-based database that contains 45120 peptides with annotated activities retrieved from multiple sources, from which 2004 are hemolytic peptides [16]. In recent years, some efforts have been made to utilize the information from these databases and predict the hemolytic activity of peptides using machine learning (ML) algorithms [1,2,4,8,9]. However, to our knowledge, no effort has been made to explore the feature space of hemolytic peptides using network science to elucidate the defining characteristics that make certain potential therapeutic peptides hemolytic.

Network science has been previously applied to successfully model many real-world systems [17]. For instance, the ‘small world-model’ introduced by Watts and Strogatz [18] has helped to understand the way different areas of the brain communicate with each other [19]; and more recently, during the Covid-19 pandemic, network science concepts were applied to develop strategies to lower the spread of infection [20,21]. Concerning therapeutic peptides, network science has been recently used to explore the chemical space and build prediction models for tumor-homing peptides [22] and antiparasitic peptides [23], having promising results.

Hence, following the same approach, this report aims to get insight into the chemical space of hemolytic peptides from the StarPepDB using network science and visual (interactive) data mining. Useful information can be retrieved from this strategy, including identifying/delineating the structural diversity among hemolytic peptides, most central and atypical peptides (singletons), the relationship between hemolysis and certain therapeutic activities, and identifying motifs related to hemolytic activity which can be useful when designing therapeutic peptide drugs. Moreover, relatively small subsets of hemolytic peptides can be extracted for further studies. These subsets (called “scaffolds”) have the advantage of representing the whole chemical space of hemolytic peptides but just using a fraction of the nodes of the complex network of hemolytic peptides [24].

Here, we describe for the first time the use of Half-Space Proximal Networks (HSPNs) to represent the chemical space of hemolytic peptides; such networks have been previously used only to explore the antiparasitic peptide space [23]. These networks possess many advantages as they generate highly connected but sparse networks that contain the minimum spanning tree as a sub-graph [25,26]. Moreover, these networks do not strictly need a pairwise similarity threshold (t) between peptides for the construction of informative networks, as is the case of Chemical Space Networks (CSNs) described in other studies. Nevertheless, despite a cutoff value t is not mandatory for HSPNs, it might affect the representativeness of the scaffolds. Hence, we compared HSPNs without a cutoff value (namely t = 0.00) with networks constructed using the same parameters but generated with their optimal similarity cutoff value. Other comparative analyses were also conducted to study the HSPN construction and visualization phase involving the use of distance metrics and centrality measures, respectively, while for extracting a representative subset of hemolytic peptides from the HSPNs, different centrality measures and global and local alignments were also evaluated.

Although, previous studies based on network science have employed the Euclidean distance as the default similarity measure metric; it is suggested that the use of different (dis)similarity measures allows the codification of orthogonal information. Hence it should not be assumed that only one measure is the best suited for calculating the similarity between objects, especially in high dimensional space [27,28,29]. For this reason, we evaluated five different two-way (dis)similarity measures for constructing such HSPNs: (1) Angular Separation, (2) Bhattacharyya, (3) Chebyshev, (4) Euclidean and (5) Soergel. Finally, by using community information from these networks and an alignment-free method for motif discovery, we reported new putative motifs that hallmark hemolytic peptides along with their further enrichment on external datasets to validate their significance.

2. Materials and Methods

The overall workflow consists of four stages: (i) Metadata network visual mining, (ii) HSPNs generation and analysis, (iii) scaffold extraction and exploration, (iv) motif discovery and enrichment (Figure 1). The first step (Section 2.2.1) involves the generation of metadata networks (METNs) and exploration of critical features related to hemolytic peptides. The second step (section 2.2.2) consists in building HSPNs that represent the chemical space of hemolytic peptides retrieved from StarPepDB. Then the best HSPN candidates were selected based on global network descriptors for further analysis (section 2.2.4). In the third step (section 2.3), representative subsets (scaffolds) from the best HSPN candidates, built up with the optimal t value, and from their respective networks with cutoff t = 0.00 were extracted by using sequence alignment and centrality information from each peptide in the graph. Finally, the last step (section 2.5) consists in proposing new putative hemolytic motifs by using an alignment-free approach and by comparing them with reported hemolytic motifs using benchmark datasets (enrichment analysis) to further select the most representative ones. All the steps of this section were performed using the StarPep toolbox, aided with in-house python scripts and the SeqKit toolkit [30].

2.1 Datasets

The datasets used in this study to construct METNs, HSPNs, generate hemolytic motifs, and subsequently conduct motif enrichment analyses are described in Supplementary Material (file SM1). The usage of each dataset in this report is detailed below.

• StarPepDB. It is a graph database embedded in the StarPep toolbox that consists of 45120 peptides with annotated activities retrieved from 40 bioactive databases and other sources [16]. A sub-dataset consisting of 2004 hemolytic peptides was extracted from this database to generate HSPNs, METNs and discover new hemolytic motifs. In addition, the complete StarPepDB was also used in the motif enrichment process to help find the most representative hemolytic motifs.
• HemoPI-1. It encompasses 552 experimentally validated highly hemolytic peptides (positive) and 552 random peptides extracted from Swiss-Prot (negative) [8]. This dataset was only used in motif enrichment analysis.
• Big-Hemo. It is a non-redundant combination of several datasets that contain either hemolytic or highly hemolytic peptides as positive samples and non-hemolytic or low hemolytic peptides as negative samples. The datasets used to generate the Big-Hemo dataset are HemoPI-2 Main and Validation [8], HemoPI-3 Main and Validation [8], HAPPENN [1], HLPred-Fuse Layer 2 Training and Independent datasets [4] and HemoNet [9]. To construct Big-Hemo, only positive samples labeled as “highly hemolytic” were retrieved from these datasets to handle the problem of lack of agreement and standardization at considering when a peptide is hemolytic or not, and the way of measuring this property, respectively [1,31]. Although HAPPENN dataset contains positive samples not labelled as highly hemolytic, its positive samples were also included in Big-Hemo in order to gain more diversity and a better representation of hemolytic peptides. Thus, this dataset was addressed to evaluate whether our novel motifs are enriched in highly hemolytic peptides, which are more concerning when designing therapeutic peptides. In addition to redundancy removal, peptides containing ‘X’ several times in a sequence and Nphe or Nleu in their sequences were also discarded. The resulting Big-Hemo dataset contains 2196 highly hemolytic peptides. Like HemoPI-1 dataset, Big-Hemo was also used for motif enrichment analysis.

Figure 1. Workflow overview of the experimental procedure. Figure created with Inkscape [32].

Figure 1. Workflow overview of the experimental procedure. Figure created with Inkscape [32].

2.2 Network Generation and Analysis

2.2.1. Metadata Networks (METNs)

A Metadata Network (METNs) is an unweighted pseudo-bipartite graph defined as $F=(V,E)$, where $E\left(F\right)$ is the set of edges of the graph and $V\left(F\right)$ is the set of nodes or vertices which comprises two classes: hemolytic peptides and metadata information (e.g., origin and function of peptides). In these networks, peptide nodes are adjacent to their corresponding metadata nodes. For instance, if a peptide is hemolytic, an edge will connect this peptide node to the hemolytic metadata node. However, METNs are not fully bipartite graphs since in the last ones the nodes belonging to the same class cannot be adjacent [33], whereas METNs can set edges within the metadata class as long as one node is hierarchically related to another. For instance, for the “Function” metadata, “Toxic to mammals” is hierarchically connected to “Hemolytic”, thus an edge connects these two metadata nodes.

StarPep toolbox allows the easy construction of METNs, which help to get insight into the related data associated with the hemolytic peptides. A Database METN, for instance, shows the databases where each hemolytic peptide has been reported by connecting each peptide node to its corresponding database nodes. This information is useful to get an overview of the most populated databases with hemolytic peptides, to analyze peptide redundancy in different databases and also to detect what peptides are uniquely reported in particular dataset, etc. Hence, we created four METNs based on different metadata information: database, function, origin and target. The peptide class of $V\left(F\right)$ was the set of 2004 hemolytic peptides from StarPepDB.

2.2.2. Half-Space Proximal Networks (HSPNs)

HSPNs, are weighted graphs defined as $G=(V,E)$ where $V$(G) represents the set of nodes (hemolytic peptides) and $E\left(G\right)$ represents the set of edges. The nodes are characterized by vectors whose components are values of sequence-based molecular descriptors (MDs), whereas the edges link nodes in a pairwise manner following the subsequent steps:

• A (dis)similarity measure is calculated for each pair of nodes using the vectors of peptide features. Then these values are normalized (min-max normalization). This forms a symmetric similarity matrix $M$ of size $n\times n$ where $n$ represents the number of hemolytic peptides and $M_{i,j}$ represents the similarity score between the nodes $V_i\left(G\right)$ and $V_j\left(G\right)$, being 1 the highest similarity value and 0 the lowest. Then a rule called Half-Space Proximal (HSP) test [26] is applied to construct the HSPN, which is a strongly connected but sparse network [25], that preserves the number of nodes while containing a relatively low number of edges compared to the counterparts, CSNs [25].
• Finally, a threshold or cutoff value $t$ can be applied to the weighted edges to further reduce the density of the graph by removing edges whose similarity values are lower than $t$. This helps to study the topology of the resulting graphs and subsequently find the best representative network of the chemical space occupied by hemolytic peptides. It is worth mentioning that for the construction of HSPNs, using a $t$ value is not mandatory.

HSPNs were constructed as follows. From the 45120 peptides found in StarPepDB [16], 2004 peptides with known hemolytic activity were retrieved using the query option of StarPep toolbox [25]. Redundancy in the peptide sequences was removed using Smith-Waterman local alignment [34] and BLOSUM-62 substitution matrix [35] considering at least 98% sequence identity, resulting in 1647 peptides (SM1.1.3). Then MDs were calculated for each peptide sequence and an unsupervised feature selection was performed, removing near constant peptide features using Shannon entropy (threshold 10%), whereas redundant features were removed using Spearman correlation coefficient (threshold 0.8%). Then all the remaining peptide features were selected for generating the networks. See reference [25] for a detailed description of the peptide feature extraction method.

Regarding the (dis)similarity measures, HSPNs were constructed using Angular Separation (AS), Bhattacharyya (Bh), Chebyshev (Ch), Euclidean (Eu), and Soergel (So) measures. Their formulae and properties are stated in Table 1. We tested several measures since previous studies demonstrated that different distance measures can codify orthogonal information; thus, not necessarily the Euclidean distance might be the best suited for a specific application [27,28].

In addition, to explore the behavior of HSPNs when varying the value of $t$, 11 different cutoffs were applied for each metric: 0.00 and from 0.50 to 0.95 in steps of 0.50, resulting in a total of 55 HSPNs available at SM3 (i.e., 11 networks for each metric). We applied these cutoffs, since a previous study showed that when constructing HSPNs, most of the global parameters barely changed when the similarity cutoff t ranges between 0.00 and 0.45 [23].

Finally, since several combinations will be generated in the following steps, we will use the following notation when referring to a specific network: “cutoff (t)_metric”. For instance, for a network generated with a t = 0.00 using the metric Angular Separation, its corresponding name will be: 0.00_AS.

2.2.3. Network Visualization

For METNs, Betweenness Centrality [36] was calculated and the size of metadata nodes was proportionally projected according to the corresponding centrality value. Database, Function, Origin and Target METNs were visualized by coloring their metadata nodes: aquamarine, yellow, light violet and green, respectively. Metadata nodes related to hemolytic activity (i.e., toxic, toxic to mammals, hemolytic, Red Blood Cells) were colored red. On the other hand, all peptide nodes had the same size and were colored blue-green for all the METNs. Database and Function METNs displayed their most central metadata nodes numbered. Finally, Force Atlas 2 layout algorithm was always used to visualize METNs [37].

For HSPNs, the nodes were clustered using the Louvain method [38], and the Hub-Bridge centrality (HB) measure was calculated for each node. Finally, to better visualize the networks, we colored the nodes according to the cluster they belong to, and the node size was set to be proportional to its HB centrality value using the Bezier interpolator. Finally, we applied the Fruchterman Reingold layout algorithm [39]. The resulting METNs and HSPNs were exported as GraphML files and further visualized with Gephi 0.9.7 [40].

2.2.4. Selection of the best HSPNs

Using Gephi 0.9.7, the following global network parameters were retrieved for each HSPN: number of edges, modularity, density, average clustering coefficient (ACC), number of clusters/communities, singletons GC (nodes disconnected from the giant component), singletons D0 (nodes of degree zero), diameter, average path length, average degree and the probability of $k$ (degree distribution). These features were used to study the behavior of the networks and select the best representations for each (dis)similarity measure (five networks in total). The best networks with their optimal cutoff value $t$ were then used for scaffold extraction.

2.3 HSPNs Scaffold Extraction and Analysis

This step aims to retrieve representative subsets of the hemolytic peptide space. The five best networks obtained in Section 2.2.4 and the five networks with $t=0.00$ (one for each similarity measure) were used to build the scaffolds. The following steps were applied for each of these networks:

The selected HSPNs were generated again following the steps of Section 2.2.1, but now only the corresponding cutoff value $t$ for each network was applied. We also calculated the Harmonic centrality (HC) for each node. After that, we applied the scaffold extraction method (integrated into the StarPep toolbox), which retrieves the most central and unique hemolytic peptides by ranking each peptide in decreasing order regarding their centrality and then redundant sequences were removed as follows: if a pair of sequences have a percentage identity higher than a certain cutoff value $s$1, the least central peptide of the pair will be removed. Finally, the resulting scaffolds of peptide sequences were exported as fasta files. This method generally assures extraction of the most representative peptides from all the centrality ranges but removes sequence redundancy.

Following the same notation for HSPNs, for naming the scaffolds, we inherit the name of the parent network followed by the centrality measure, alignment type, and cutoff $s$ value: “cutoff t_metric_centrality_alignment_cutoff s”. For instance, a scaffold extracted from the network 0.00_AS using harmonic centrality, local alignment and a cutoff s = 0.80 would be named as: 0.00_AS_HC_L_0.80.

In this experiment, we varied the type of centrality measure, the sequence alignment type and the cutoff value of percentage identity s. We used Hub-Bridge (HB) or Harmonic centrality (HC), and Needleman-Wunsch global alignment (G) [41] or Smith-Waterman local alignment (L) [34] were used for sequence comparison, both with BLOSUM-62 substitution matrix. Moreover, we tested various cutoff values $s$ ranging from 0.40 to 0.90 in steps of 0.10. As a result, for each of the ten selected networks, we generated 24 different scaffolds using the combinations described above. In total 240 scaffolds were obtained (see SM4).

Scaffold Comparison by Metric. For all the scaffolds generated from networks where the cutoff $t=0.00$, the Jaccard similarity coefficient (JSC) was calculated between scaffold pairs created with the same parameters but differing in their metric. For instance, 0.00_AS_HB_G_0.40 vs. 0.00_Bh_HB_G_0.40. JSC is defined as the number of elements of the intersection of sets A and B divided by the number of elements of the union of those sets [42]. We calculated this distance to assess the similarity between scaffold pairs generated with different parameters.

Scaffold Comparison by Cutoff t, Type of Alignment and Centrality Measure. Using scaffolds constructed from orthogonal metrics (section 2.4.1), we compared the effect of varying the $t$ value when generating networks. For this task we calculated the JSC between scaffold pairs created with the same parameters but differing in their $t$ value. For instance, 0.00_AS_HB_G_0.40 vs. 0.90_AS_HB_G_0.40. Similarly, the same approach was followed to evaluate the effect of the type of alignment and centrality measure in the representativity of the scaffolds.

All pairwise scaffold comparisons were conducted using SeqKit toolkit to extract peptide IDs and then the relation between sets was obtained using https://bioinformatics.psb.ugent.be/webtools/Venn/.

2.5. Motif Discovery and Enrichment

2.5.1. Motif Discovery

Motif discovery was performed using the alignment-free method STREME (short for Sensitive, Thorough, Rapid Enrichment Motif Elicitation), which finds ungapped motifs enriched in input sequences compared to control sequences providing a statistical significance for each motif [43]. To generate a diversity of potential new hemolytic motifs, we employed the community information of HSPNs with cutoff $t=0.00$ generated with the metrics: Angular Separation, Chebyshev and Euclidean. The following steps were performed for each of the networks:

• Using the StarPep toolbox we extracted the sequences of peptides belonging to each cluster (community) and saved them as fasta files. Then these files were used as input sequences for motif discovery. For control sequences, we let the method use shuffled input sequences. Since our peptides contain non-standard AAs, we provided a customized alphabet (SM5.1.1). Motifs ranging from 3 to 6 letters, at least 20% present in the input sequences and with a p-value lower than 0.05 were retrieved.

Similarly, the same steps were applied to retrieve motifs from the file containing 1647 non-redundant hemolytic peptides from StarPepDB. Then, motifs resulting from this process were combined with hemolytic motifs reported in the literature [1,8], and duplicated motifs were removed (see SM5.1.2).

2.5.2. Motif Enrichment

Motif Enrichment was conducted using SEA (Simple Enrichment Analysis) from the MEME suite (https://meme-suite.org/meme/tools/sea). Hemolytic motifs found in section 2.5.1 and motifs reported in the literature [1,8] were employed to assess whether they are enriched in benchmark databases. These motifs were evaluated in HemoPI-1, StarPepDB and Big-Hemo datasets. Sequences labeled as non-hemolytic were used as control sequences in enrichment analysis on HemoPI-1; sequences not having “hemolytic” metadata were used as control when the StarPepDB was used. For enrichment analysis in the Big-Hemo dataset, input sequences were shuffled by the SEA algorithm and used as control sequences. In addition, sequences with length less than three AAs were discarded (SM1.1.3). Finally, those motifs that are statistically significant in all three datasets were kept.

3. Results and Discussion

3.1. Metadata Networks (METNs)

Database METN. Most hemolytic peptides of the StarPepDB come from the SATPdb [12], Hemolytik [14], DBAASP [15], UniProt [44], DRAMP [45] and CyBase [46] databases that are the six most central nodes in Figure 2A. The majority of peptides are shared by SATPdb, Hemolytik, DBAASP and DRAMP, whereas CyBase contains more unique sets of peptides. It might be because CyBase mainly focuses on collecting information about specific types of proteins, cyclic proteins which have shown to possess important advantages such as higher stability and binding affinity compared with linear peptides [47].

In addition, SATPdb has the highest betweenness centrality and node degree value since it is connected to 1817 hemolytic peptides. On the contrary, the databases having the least number of hemolytic peptides are NeuroPep [48], Defensins [49] and Bagel 2 [50] which have node degrees of 4, 2 and 1, respectively. Overall, the Database METN can be helpful when searching for the most important databases regarding peptide hemolytic activity as well as the most unique and most specialized databases.

Function METN. When designing therapeutic drugs, understanding other activities associated with hemolytic peptides can be a good starting point for inferring possible mechanisms of action or chemical characteristics of peptides that might be related not only to certain therapeutic activity but also with hemolysis. A Function METN can be a fast and easy approach to tackle this question by using the StarPep toolbox. Figure 2B shows a Function METN of the 2004 hemolytic peptides reported in the StarPepDB. Evidently, the most central activities are “hemolytic”, “toxic” and “toxic to mammals” since the peptides of study are hemolytic and the aforementioned metadata nodes are hierarchically related (colored red in Figure 2B with centrality ranks: 1, 3 and 5, respectively). However, most of these peptides are also related to antimicrobial activity and hierarchically related metadata: antibacterial, anti-Gram positive, anti-Gram negative, antifungal, etc. In fact, these metadata comprise the nine most central nodes in the Function METN.

Since the main target of antimicrobial peptides (AMPs) is the bacterial cell membrane which is disrupted by several reported modes of action [31], it might be feasible that similar modes of action can also target and disrupt human cells, specifically RBCs. Many studies have proposed that due to the positive charge of many AMPs, they can selectively disrupt negatively charged membranes of bacteria while not affecting the neutral membranes of mammals [51,52]. However, it has been demonstrated that several AMPs (some with high antimicrobial activity) can also disrupt mammalian cells as well, causing hemolysis in RBCs [31,53]. In fact, Function METN shows that 94.46% of the 2004 peptides that comprise the hemolytic space, have both antimicrobial and hemolytic activity.

Origin METN. This type of METN helps to easily identify the origin of hemolytic peptides, whether they are synthetic or isolated from living organisms. Figure 3A shows the complete Origin METN in the dashed box. The central part of the METN was zoomed in and depicted in the center of Figure 3A. Looking at the complete Origin METN three distinctive regions can be observed, an outer ring, a middle ring and a central network. The outer ring represents peptides isolated from living organisms but have not been chemically synthesized. For instance, the peptide StarPep_06954 [54] whose metadata origin node corresponds to only Caenorhabditis elegans. The middle ring represents peptides with nodes of degree zero. This means, peptides that do not have an origin label in the StarPepDB.

On the other hand, the central network shows peptides that have only synthetic origin (the most central blue-green nodes) and peptides isolated from living organisms that have also been chemically synthesized (nodes connected to the central violet metadata node and also connected to radial violet nodes). Radial violet nodes connected in a chain-like way represent hierarchical taxonomic ranks that are related to species from which a particular peptide was obtained. For instance, the subsequent metadata nodes are connected in the following manner Urochordata->Ascidiacea->Pleurogona->Stolidobranchia->Pyuridae->Halocynthia->Halocynthia aurantium. The H. aurantium metadata node is then connected to 6 peptide nodes isolated from that species.

Over half of the hemolytic peptides (1060) are of synthetic construct, whereas the rest are isolated from various organisms. From the top 20 most central origin metadata nodes (synthetic construct not included), half of them belong to the class Amphibia. This is expected because most of the hemolytic peptides in the StarPepDB are antimicrobial (Figure 2B) and a significant part of them have been isolated from frogs and toads since it has been known that they can produce broad-spectrum AMPs in their granular glands in the skin as a defense strategy [55,56,57].

Target METN. An outer ring and a central network can be observed in this METN (Figure 3B). The outer ring of peptides seen in the dashed box are peptides that do not have a metadata node related to a target. This metadata network works in the same fashion as the Origin METN, where chain-like nodes represent the hierarchical taxonomic ranks, but instead of representing the origin of the peptide, it displays the target of the peptide i.e., the species/cell type in which a certain peptide activity has been evaluated. Evidently, the main target is human erythrocytes (colored red in Figure 3B) since we are exploring the hemolytic peptide space. Other central targets include Escherichia coli, Staphylococcus aureus, Pseudomonas aeruginosa, Bacillus subtilis and Candida albicans. They are among the six most central metadata nodes in this METN. It shows that several of the hemolytic peptides have been evaluated as potential AMPs in important human pathogens such as P. aeruginosa which has become a real concern in hospital-acquired infections due to drug-resistance appearance [58].

GraphML files of METNs and the descriptor information from each node are available at SM2.

3.2. Half-Space Proximal Networks (HSPNs)

The properties of the HSPNs were studied based on their global network parameters consisting of the number of edges, modularity, density, average clustering coefficient (ACC), number of communities and singletons, among others. Such statistics can provide a good picture of the topology of the graphs and help selecting networks with the cutoff t that better projects the chemical space of hemolytic peptides.

Our results are consistent with another study that showed that there was little change in the global network parameters when networks are created within the cutoff t range 0.00 – 0.45 [23]. This is because of the highly low number of edges that are removed within this range. In fact, on average, the number of removed edges at $t=0.50$ correspond to the 1.9% of the initial edges when $t=0.00$ (See SM3.6).

Moreover, it can be observed that networks generated by different metric measures address differently the similarity between peptides (Figure 4). Based on their behavior, the networks used in this study can be roughly grouped into three classes: Class I: Angular Separation; Class II: Bhattacharyya, Euclidean and Soergel; and Class III: Chebyshev. The influence of the metric measure in the global parameters of the networks is provided below. All global network parameters calculated for each metric are provided in SM3.

Modularity. This is a measure of the network connectivity which indirectly represents how well-defined communities are in the graph and is associated with the number of communities. Graphs generated with Angular Separation (AS) initially possess higher modularity values compared to the other metrics; however, the modularity keeps relatively low at higher t values (0.550 at t = 0.95) whereas the other four metrics increase their modularity to values near 1. On the other hand, Chebyshev (Ch) networks show the lowest modularity at low cutoff values, but then it increases to high values comparable with Soergel (So), Euclidean (Eu) and Bhattacharyya (Bh). So- and Eu-derived networks have a quite similar behavior in the whole range of t values, whereas Bh networks initially behave similar to Eu and So networks, but then diverge at t = 0.70 (Figure 4A). An adequate selection of modularity is important since highly sparse networks with an elevated number of communities would not provide useful information as several resulting communities would be just artifacts.

Density. It shows the ratio between the edges present in the network and the maximum number of possible edges. Similarity networks have been shown to have an inversely proportional relationship between similarity threshold (t) and density [25,60]. The same pattern is observed for all metrics, but with some notable variations. Here, we can identify three behaviors according to the three class of metrics. AS networks have the highest density in the entire range of t, whereas Class II metrics (i.e., Bh, Eu, and So) have the lowest density until t = 0.70. On the other hand, Ch networks not only have an intermediate initial density but also show the biggest variation of density along the whole range of t (Figure 4C). In order to select adequate networks, we should choose graphs that are neither too dense nor too sparce since the former would hamper retrieval of useful information whereas the latter would lose information [61]. Density values below 0.20 are desired as they allow to properly understand the network while preserving high modularity. Particularly, HSPNs are suited because they have the intrinsic characteristic of showing low densities. In fact, the highest density value in this study corresponds to 0.020 (0.00_AS network).

Average Clustering Coefficient (ACC). This measures the connectivity of the network and it has been previously studied on molecular similarity networks varying the similarity cutoff t. One study showed that the ACC maximum peak correlates with the best clustering outcome and is a good indicator for finding the appropriate value of t [60]. In our study, three behaviors related to the metric class can be observed again. AS networks have the highest ACCs in the whole range of t with its local maximum at t = 0.95. On the other hand, Ch networks start with very low ACCs and get increased at t = 0.65 reaching its maximum peak at t = 0.90. Finally, Class II metrics have the lowest ACCs in the entire range of t with their maximum peaks at 0.70 (Eu), 0.80 (So) and 0.85 (Bh) (Figure 4E).

Communities and singletons. The number of communities determined with the Louvain method, the number of singletons D0 (nodes of degree zero) and the number of singletons GC (nodes disconnected from the giant component) were calculated to select the networks with the most reasonable values of these parameters. When t = 0.00, HSPNs have the minimum spanning tree as a subgraph, this implies that at this t value all nodes are connected. In other words, no singletons D0 nor singletons GC are found. Regarding the number of communities at t = 0.00, all metric networks showed similar values (on average 8 communities). At higher t values, the number of communities and singletons D0 increase dramatically for all the metric networks, except for AS networks (Figure 4B–D). This is expected as more edges are removed, more nodes are isolated, and now singletons are counted within the communities. Hence, an appropriate t value should be selected that comprises an equilibrium between singletons (atypical peptides) and communities that reflect a real chemical relationship.

Other global network parameters were also calculated to characterize the networks, such as the diameter of the graph (Figure 4F), the average path length (APL) and average degree (See SM3.6). In order to find the best t value for each metric network, we should look for a compromise between the best parameter value for each descriptor i.e., networks with low density, with neither too many clusters (< 20) nor too many singletons (~15-30), retaining high ACC and high modularity. The global descriptors of the selected networks with their best cutoff value t and their respective networks constructed with t = 0.00 (10 networks in total) are shown in Table 2.

Finally, we calculated the probability of k (also known as the degree distribution) for each of the selected networks (Figure 5). Overall, all networks show a right-skewed bell-shaped distribution with high probability of intermediate node degrees. Evidently, plots on the left (t = 00) show a probability of zero for singletons (k = 0) whereas plots on the right (best value t) tend to have a higher probability when k = 0. In addition, plots with the best t value have smaller maximum degrees (as well as the average degree) compared with same-metric networks at t = 0.00. Thus, when comparing networks with the same metric but varying the cutoff value (t = 0.00 vs. best cutoff t), it seems both retain a similar degree distribution. However, when comparing networks with different metrics we can get marked differences. AS networks tend to have a wider distribution range and a higher average degree whereas Ch networks show intermediate values, and networks constructed with Class II metrics show a similar distribution shape among them and have the lowest distribution ranges and average degrees of all metrics. Figure 6 shows the graphical representation of the 10 selected HSPNs.

3.3. HSPNs Scaffolds

A total of 240 scaffolds were extracted from the 10 HSPNs selected in section 3.2 (SM4.1). To better understand the effect of the centrality measure, type of alignment and cutoff value s when constructing the scaffolds, several pairwise similarity comparisons between scaffolds were carried out using the JSC between sets of hemolytic peptides [42].

Metric Comparison. We compared the type of metric measure used to build the parental networks of the scaffolds. For this comparison, scaffolds (t = 0.00) built with the same combinations of centrality, alignment, and cutoff s but with different metrics were evaluated (SM4.3.1). Each pair of scaffolds is represented as a point in Figure 7. In all plots of Figure 7 when s ≥ 0.60, all scaffold pairs constructed with Class II metrics (i.e., Bh, Eu, So) show the highest similarity percentage compared with the pairs from other combination of metrics. Moreover, scaffold pairs in which one of them is extracted by the AS metric show the smallest similarity percentage at almost any cutoff value s. On the contrary, scaffolds selected with Ch metric have an intermediate similarity percentage when compared with scaffolds extracted by other metrics. These results agree with section 3.2 which showed that the five metrics tend to have three types of behavior (three classes of metrics). The density (Figure 4C) and the degree distribution (Figure 5) of the networks with different metrics are the global descriptors most correlated with the results from the percentage similarity among scaffolds. Thus, it is possible to reduce the number of highly-similar scaffolds by using only those HSPNs with the metrics that mostly differ in the aforementioned global network parameters. In this case, Class II metrics: Bh, Eu, and So are the metric measures with the most similar behavior since they produce similar networks and scaffolds. Therefore, it was decided to conduct the following analyses using only one of the metrics of Class II: Euclidean. This metric was chosen since it is the default metric used in other studies [22,23], and it would be advantageous to compare its performance with the other metrics not previously used in this type of study. Overall, this step allowed us to reduce the redundancy in the scaffold representativity from 240 to 144 scaffolds (SM4.2).

Cutoff Comparison. A cutoff value t is not mandatory when constructing HSPNs since at t = 0.00, these networks already have low densities under 0.20. However, the topology, characterized by global network features, tends to vary when varying t as was demonstrated in section 3.2. Thus, it is important to evaluate the effect of selecting a cutoff value (or not) when constructing representative scaffolds of the chemical space. The JSC was calculated between pairs of scaffolds extracted by using the same metric but at different cutoff values (t = 0.00 vs. best t value), see SM4.3.2 (Figure 8).

A marked difference was observed when these scaffold pairs were constructed with different types of centralities. Scaffolds constructed with HB centrality (Figure 8A,B) tend to have more unique peptides at low s values and the number of common peptides between scaffold pairs tend to increase when s increases. A similar pattern was observed in Figure 7. However, when the same scaffolds are constructed replacing HB centrality with HC centrality all scaffold pairs tend to share more than 89.50% of peptides regardless of the value of s (SM4.3.2.2) (Figure 8C,D). Furthermore, the same patterns are preserved when any alignment type is applied. Hence, when generating scaffolds using HC centrality, it is unnecessary to first find the best t value for the parental networks since similar scaffolds will be obtained using networks with t = 0.00.

Alignment Comparison. A clear pattern can be observed when extracting scaffold either using global or local alignment (Figure 7 and Figure 8). In general, local alignment tends to discriminate more strongly at low s values than global alignments. Hence, scaffold pairs extracted with local alignment at such low s values have a lower similarity percentage than the analogue scaffold pairs extracted using global alignment.

In addition, when comparing the similarity percentage of scaffold pairs extracted using the same parameters but differing the alignment type, the same behavior was observed independently of the metric, type of centrality or the t value used, see Figure 9. Scaffold pairs differing only in their alignment type tend to have a low percentage of similarity at low s values, which might indicate that these methods capture the similarity between peptides differently. However, when analyzing the proportion of unique peptides between these scaffold pairs, it is clear that scaffolds extracted using local alignment are practically a subset of scaffolds extracted when using global alignment. In fact, the average number of unique sequences in local scaffolds when comparing them with their global counterparts at any cutoff s is 16.19 (SM4.3.3). An example is provided for the scaffold pairs: 0.00_AS_HB_G_0.40 and 0.00_AS_HB_L_0.40 (Figure 10).

Centrality Comparison. Pairwise comparisons of the scaffolds constructed using the same parameter but changing the centrality measure show a trend like the pairwise comparisons presented before (SM4.3.4). This implies that the type of centrality used to extract the scaffold will affect in the sequences that are removed/retained, especially at low s values.

On the other hand, when comparing centrality measures, JSC between scaffold pairs extracted from networks with best t value tend to be higher than JSC from scaffold pairs from networks with t = 0.00. This pattern is clearer at low s values (Figure 11).

All scaffolds presented in this section can be used in many applications. For instance, they can be used as training datasets for both ML-based and Multi-Query Similarity Searching (MQSS) prediction models of hemolytic peptides. The advantage of using these datasets is that they store information of central and important peptides as well as outliers or atypical hemolytic peptides while avoiding overrepresentation of certain peptide classes (sampling bias). Each scaffold bears a unique type and amount of information of the hemolytic peptide space and one scaffold can be more suitable than another depending on the scaffold’s use. Scaffolds extracted at low cutoff s values tend to cover fewer peptides of the original space, whereas higher s values capture more information of the space but peptide overrepresentation might be present. Even scaffolds with s = 0.90 have on average 71.48% coverage of the total hemolytic peptide space, which is still an adequate equilibrium between coverage and amount of information. Figure 12 depicts an example of the scaffold coverage when varying the cutoff s.

3.4. Hemolytic Motif Discovery and Enrichment

Motif Discovery. Peptides from each community were used as input sequences to uncover new hemolytic motifs within the communities’ diversity by means of STREME, an alignment-free method [43]. Table 3 shows a sample of the 42 new motifs discovered using clusters of HSPNs (t = 0.00) created with different metrics. 12 motifs were found from 6 clusters of the network 0.00_AS, 14 motifs were discovered from 4 clusters of the network 0.00_Ch and 16 motifs from 5 clusters were discovered using the network 0.00_Eu. The three metrics commonly detected only 4 motifs: GLP, MFTKL, ERBADE and VCTRN. It is worth mentioning that several other motifs were similar but not identical such as: GLP/GLPV or VGGTCN/GGTCN. In addition, 15 motifs were discovered not considering the community diversity by using all 1647 hemolytic peptides as input sequences. All these motifs were grouped as HSPNs motifs. After removing duplicated motifs, 50 HSPNs motifs were discovered (SM5.1.2).

Two previous reports on ML models for predicting hemolytic activity of peptides have also reported hemolytic motifs, namely: HemoPI [8] and HAPPENN [1]. HemoPI reported 21 motifs extracted using MERCI software that were enriched in positive sequences from HemoPI-1 and HemoPI-2 datasets, whereas HAPPENN motifs resulted by looking for the 20-top motifs found exclusively in the positive dataset of HAPPENN. No HSPN-derived motifs were found among the reported ones. To generate a unique list of non-redundant hemolytic motifs, HSPNs motifs were combined with the previously reported ones resulting in 91 putative motifs. Then similar motifs were combined into consensus motifs resulting in 57 non-redundant motifs (SM5.1.2 and SM5.1.3).

Motif Enrichment. To identify and validate the most representative hemolytic motifs and remove some artifacts from the 57 potential hemolytic motifs, we conducted enrichment analyses using SEA method on three different datasets: HemoPI-1, StarPepDB and Big-Hemo (SM5.2). Motifs not reported as significant in at least one dataset were removed. The resulting 47 hemolytic motifs sorted by the average enrichment ratio of all datasets are presented below (newly discovered motifs by HSPNs are shown in red):

• MFTLK, ALKAIS, GTCN, WKSFJK, VCGETC, WKK, AKKAL, GETCV, CYCR, LKKL, CVCV, ISWIK, RFC, LHTA[KL], FLHSAK, CSW, LWKT, FLGTI, GAVLKV,PGC, KKILG, KITK, KHI, LGKL, KWK, VNWK, K[GT]AGK, VCT, ALW, SWP, HIF,LLKK, [VI]LDTJ, CRR, KLL, JGKL, FKK, GAIA, VLK, GLP, PKIF, GKEV, GTIS, AAAK, GCS, IAS, MAL (Table 4).

These motifs might be involved in the mechanisms of action of hemolytic peptides as well as antimicrobial activity but further studies are needed to corroborate this assumption. Another possible use of these motifs can be as a toxic signature, where proteins containing some of these motifs could be attributed to a relatively high hemolytic activity in comparison with proteins with few or none hemolytic motifs. Table 5 shows an example of three pairs of peptides whose hemolytic activity is related to the number of hemolytic motifs present in their sequences.

4. Conclusions

Positive endpoints are commonly evaluated in peptides to better understand a specific therapeutic activity; however, getting insight into negative endpoints such as hemolysis should be equally important. In this study, the exploration of the chemical space of hemolytic peptides through a synergic combination of network science and interactive data mining resulted in an easy and feasible way to get more insight into the features that characterize this type of peptides. For instance, the hemolytic activity was almost exclusively associated with AMPs, from which the majority are of synthetic construct or isolated from frogs and toads. Moreover, scaffolds and motifs extracted from the hemolytic peptide space can aid future studies related to description of hemolytic peptide families, assessment of the mechanisms of hemolysis, and training ML- and MQSS-based prediction models. The main asset of this approach is that it can easily be extrapolated to many other peptide types such as AMPs, industrial enzymes, therapeutic hormones, etc.