Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 journal homepage: www.elsevier.com/locate/csbj HPIPred: Host–pathogen interactome prediction with phenotypic scoring Javier Macho Rendón, Rocio Rebollido-Ríos 1, Marc Torrent Burgas ⇑ Systems Biology of Infection Lab, Department of Biochemistry and Molecular Biology, Universitat Autònoma de Barcelona, Spain a r t i c l e i n f o a b s t r a c t Article history: Protein-protein interactions (PPIs) are involved in most cellular processes. Unfortunately, current knowl- Received 15 August 2022 edge of host-pathogen interactomes is still very limited. Experimental methods used to detect PPIs have Received in revised form 9 November 2022 several limitations, including increasing complexity and economic cost in large-scale screenings. Hence, Accepted 10 November 2022 computational methods are commonly used to support experimental data, although they generally suffer Available online 21 November 2022 from high false-positive rates. To address this issue, we have created HPIPred, a host-pathogen PPI pre- diction tool based on numerical encoding of physicochemical properties. Unlike other available methods, Keywords: HPIPred integrates phenotypic data to prioritize biologically meaningful results. We used HPIPred to Protein-protein interaction Host screen the entire Homo sapiens and Pseudomonas aeruginosa PAO1 proteomes to generate a host- Pathogen pathogen interactome with 763 interactions displaying a highly connected network topology. Our predic- Bacteria tive model can be used to prioritize protein–protein interactions as potential targets for antibacterial Pseudomonas aeruginosa drug development. Available at: https://github.com/SysBioUAB/hpi_predictor. PAO1 Ó 2022 The Author(s). Published by Elsevier B.V. on behalf of Research Network of Computational and interactome Structural Biotechnology. This is an open access article under the CC BY-NC-ND license (http://creative- commons.org/licenses/by-nc-nd/4.0/). 1. Introduction rithms to predict host-pathogen PPIs is still far from optimal due to the lack of validated experimental datasets [8]. As a result, many During infection, pathogen proteins play a crucial role in re- of these algorithms display high false-positive rates [16], making wiring multiple biochemical processes in the host, ultimately them unusable for PPI discovery in the laboratory. To address these allowing the pathogen to attach to and invade host cells [1–3]. In limitations, we have developed an algorithm to predict host- return, the host uses sophisticated defense mechanisms against pathogen PPIs that uses numerical representations of proteins pathogens. Most of these processes are mediated by protein–pro- based on the physicochemical properties of their amino acids. To tein interactions (PPIs) [4]. Several detection methods, such as improve robustness, individual model predictions are then com- yeast two-hybrid, pull-down assays, or coimmunoprecipitation, bined into a consensus interactome, which is finally integrated are commonly used to identify novel PPIs [5–7], although only a with phenotypic data collected from infection-related databases, small fraction of the PPI space has been characterized so far [8]. allowing us to provide a ranked score for each interaction. In order to address the shortage of validated PPIs, in silico methods are commonly used. There are many PPI predictors available, including homology- 2. Methods based [9,10], annotation-based [11,12], structure-based [13], and deep learning methods [14,15]. However, the potential of algo- 2.1. Data collection and dataset construction Abbreviations: PPI, Protein-protein interaction; CCC, Cross-correlation coeffi- Positive dataset: Host-pathogen PPIs were obtained from cient; BC, Betweenness centrality; FPR, False positive rate; ROC, Receiver-operating PHISTO [17], a database of experimentally validated interactions. characteristic; PR, Precision-recall. A total of 9.237 inter-species PPIs between 95 different bacterial ⇑ Corresponding author. strains and Homo sapiens were used as a positive dataset. It is E-mail address: marc.torrent@uab.cat (M. Torrent Burgas). 1 worth noting that 90 % of these entries belong to Homo sapiens – Present address: University of Cologne, Faculty of Medicine, and Cologne University Hospital, Department of Internal Medicine, Center for Integrated Oncology Yersinia pestis (4.069), Homo sapiens – Bacillus anthracis (3.053), Aachen, Bonn, Cologne, Düsseldorf, CECAD Cologne Cluster of Excellence on Cellular and Homo sapiens – Francisella tularensis (1.348). Then, we applied Stress Responses in Aging-Associated Diseases, Cologne, Germany. a length cut-off to remove PPIs containing any protein shorter than https://doi.org/10.1016/j.csbj.2022.11.026 2001-0370/Ó 2022 The Author(s). Published by Elsevier B.V. on behalf of Research Network of Computational and Structural Biotechnology. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). J. Macho Rendón, R. Rebollido-Ríos and M. Torrent Burgas Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 100 amino acids or longer than 2.000 amino acids, obtaining a final windows may result in information loss, while smaller windows dataset of 7.423 PPIs. Small proteins and peptides can easily corre- would be less effective to capture the environment. late to small regions in other unrelated proteins, increasing the 2. Calculation of similarity scores by cross-correlation: To probability of false positive outcomes. By filtering entries shorter determine similarity in the physicochemical profiles of proteins, than 100 amino acids (400 proteins for the pathogen proteome we calculated the cross-correlation coefficients (CCCs) between and 700 for the host interactome) we avoid false positives while the query dataset and the positive dataset by performing one-vs- losing <0.5 % of potential interactions, on average. At the end, this all pairwise comparisons, i.e., each query protein was individually dataset represents the interactome of 3.327 human proteins tested against all proteins in the positive dataset, for both host and against 2.496 bacterial proteins. pathogen proteins (Fig. 1B). To remove low-scored pairs, all com- Synthetic negative dataset: Protein libraries containing random parisons with CCC < 0.4 were removed (Fig. 1C). Since CCCs depend sequences based on the 20-amino acid alphabet were created using on the lag parameter, we tested how many pairs were recovered by a gamma distribution to fit the observed protein-length distribu- increasing the lag interval. For intervals higher than 200, no gain tion in eukaryotic and bacterial organisms. The average protein was observed (Fig. 2) so all CCCs were calculated with a lag interval lengths used were 472 for the eukaryotic proteome and 319 for of [-200, 200]. The highest CCC, i.e., {max (CCC), lag 2 [-200,200]}, the prokaryotic proteome [18]. The synthetic proteomes contained obtained for each pairwise comparison was assigned as a measure 20.000 and 3.000 proteins, similar in size to human and bacterial of similarity between the two proteins being compared, creating a proteomes, respectively. Subsequently, the same length filtering database of similar proteins. Each entry in the database represents criterion was used to remove proteins shorter than 100 or longer a pairwise comparison and includes information on the highest than 2.000 amino acids, giving a total set of 18.669 and 2.598 pro- CCC and the length of the query protein. teins for the host and bacteria, respectively. 3. Filtering the database using a synthetic negative dataset: To Negative dataset: The Homo sapiens proteome was downloaded reduce the number of false positive predictions, we introduced a from UniProt (UP000005640) and used as host proteins. The com- filtering criterion using the synthetic dataset. We evaluated the bined proteomes of Yersinia pestis (UP000000815), Francisella CCCs for all proteins in the synthetic set against host and pathogen tularensis (UP000001174), and Bacillus anthracis (UP000000594) proteins in the positive set. As the synthetic set contains only were also downloaded from UniProt and used as pathogen pro- random-sequence proteins, none of them should be considered teins. Non-interacting pairs of proteins were generated by ran- an interacting pair. The CCCs of the synthetic dataset were then domly pairing proteins from the host and the pathogen fractions, plotted against their corresponding protein lengths. As the length discarding pairs that belonged to the positive PPI dataset. The of the proteins increased, their associated CCCs decreased linearly, length filtering criterion was applied as previously described, to a correlation that was also observed with the query datasets obtain a dataset of 7.421 entries, containing 2.734 and 2.102 differ- (Fig. 1C). Hence, we determined the slope and intercept points of ent proteins from host and pathogen, respectively. parallel linear equations of the type y = ax + b that represented Query proteome datasets: Homo sapiens (UP000005640) and such a negative linear correlation between protein length and Pseudomonas aeruginosa PAO1 (UP000002438) proteomes were CCC for the synthetic dataset, so that only 1 %, 0.1 %, 0.01 %, used to build the Homo sapiens – Pseudomonas aeruginosa interac- 0.001 %, 0.0001 or 0.00001 % of the data points fell above the equa- tome. After applying the length filtering criteria, host and bacteria tions. These linear equations were calculated and averaged over proteomes were composed of 19.192 and 1.314 proteins, the five different physicochemical parameters previously respectively. described (Fig. 1C). Afterward, entries with CCC lower than or equal to the synthetic counterpart were discarded, to obtain a fil- tered database of similar proteins (Fig. 1C). 2.2. Prediction of protein–protein interactions 4. Prediction of protein–protein interactions: To predict new The main steps involved in our prediction algorithm are: (1) the interactions, each PPI from the query dataset was inspected the fol- encoding of amino acid sequences to numerical strings, (2) the cal- lowing way: the host protein in the query PPI was searched against culation of similarity scores between the query and positive data- all the pairwise comparisons in the positive database of similar sets, (3) filtering entries by synthetic negatives, and (4) the proteins and whenever a match was found, the query interaction prediction of putative PPIs based on their similarity scores (Fig. 1). was kept. The same procedure was repeated for the pathogen protein in 1. Numerical encoding of protein sequences: Each protein the subset of positive PPIs that contain the host match. If the query sequence in all datasets (positive, negative, synthetic, and query) pathogen had a match in the subset, the interaction was consid- was transformed into a numerical string by using physicochemical ered a putative PPI (Fig. 1D). We performed this search sequen- properties of amino acids, transforming the amino acid sequences tially with all the PPIs in the query dataset to obtain the into a numerical signal (Fig. 1A). The physicochemical properties predicted interactome (Fig. 1D). To generate consensus interac- are experimentally determined values for each amino acid as tomes, individual host-pathogen interactomes were predicted for included in the AAindex database [19]. To represent the main prop- the five different physicochemical properties previously described erties that contribute to protein binding, we used five different and then combined to include any PPI that had been predicted by physicochemical indices: alpha-helix propensity (GEIM800101) at least three individual models (Fig. 3A). and beta-strand propensity (GEIM800105), to represent structure and ultimately hydrogen bonding; hydrophobicity index (ZIMJ680101) to represent hydrophobic effect; isoelectric point 2.3. Phenotypic scoring of predicted PPIs (ZIMJ680104) and electron–ion interaction potential values (COSI940101) to account for the electrostatic potential. All descrip- To add additional layers of information to the predicted interac- tors were normalized from 0 to 1 to avoid biases. After numerical tions, we compared the proteins involved against several databases encoding, a moving average algorithm with a sliding window of that contain information about the infection phenotype, namely 9 positions was used to smooth the data and represent each amino BacFITbase [20], DualSeqDB [21], and PHI-base [22] (Fig. 3B). acid’s numerical value as a measure of itself and its near environ- 1. Sequence alignment against BacFITbase: We downloaded ment. Structural elements are, on average, between 5 and 10 (beta BacFITBase v1.0 (accessed 2 September 2021), a database that con- sheet) and 3–11 (alpha helix) residues long. Hence, using larger tains information on bacterial fitness, as measured by transposon 6535 J. Macho Rendón, R. Rebollido-Ríos and M. Torrent Burgas Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 Fig. 1. PPI prediction algorithm using a single model. A) Protein sequences are numerically encoded using physicochemical descriptors and then smoothed using a sliding window approach. B) Protein profiles are then compared using cross-correlation and the calculated values are plotted against protein length. C) The dataset generated is first filtered by a coefficient cut-off (Filter 1) and later by a selected filter to control for false positive tolerance (Filter 2). D) At the end, the predicted protein pairs that survived the filtering steps are created and organized in a single interactome. 6536 J. Macho Rendón, R. Rebollido-Ríos and M. Torrent Burgas Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 Fig. 2. Optimal maximum lag value for CCC calculation. To find the optimal maximum lag value, the CCCs between the Homo sapiens query proteins and the positive dataset were calculated for different lag intervals. The plot shows the number of protein pairs with a CCC > 0.4 for each maximum lag value tested for the hydrophobicity descriptor. The number of surviving protein pairs did not increase substantially after a maximum lag of 200. Fig. 3. Model combination and calculation of ranked scores. A) Model combination. B) Determination of ranked scores for the PPIs in the consensus interactome. mutagenesis. To determine the similarity between pathogen pro- percentage of identity  40 % and an E-value < 10 that have a sig- teins in the predicted interactome and entries in BacFITBase, nificant fitness score in BacFITBase (adjusted p-value  0.05). If HPIPred performs a BLAST sequence alignment between the query multiple entries are retrieved, HPIPred assigns the average fitness pathogen protein and the entire database, keeping all hits with a score, and stores the mean standard deviation. Queries with no hits 6537 J. Macho Rendón, R. Rebollido-Ríos and M. Torrent Burgas Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 are labeled as ‘‘NA”. Finally, as BacFITbase assigns the lowest fit- F þ Eh þ Ep þ P þ BC Av S ¼ ð1Þ ness scores to the most relevant proteins, the values are normal- NM ized from 0 to 1, assigning a value of 1 to the lowest fitness where AvS is the average score, F is the fitness value, Eh and Ep score reported, and a value of 0 to the highest. are the log2 fold change in expression for host and pathogen, 2. Sequence alignment against DualSeqDB: We downloaded respectively, P is the infection phenotype, BC is the host between- DualSeqDB 1.0 (accessed 2 September 2021), a database that con- ness centrality, and NM is the number of non-missing values. Then, tains information on gene expression changes in bacterial infection a phenotypic weight (PW) weight was calculated to consider the models. Changes in gene expression are represented as the log2 number of missing values (NA) in the previous formula. Hence, a fold change, as measured by dual RNA-Seq experiments. HPIPred PW of 5 means no missing values, and 0 that no values were performs a protein sequence alignment between each query pro- reported for that specific PPI. To account for both the average score tein and DualSeqDB, for both the bacterial and host fractions, keep- and the confidence weight, we calculate a normalized ranked score ing those hits with a percentage of identity  40 % and an E- (RS) (Fig. 3B): value < 10 that have a significant expression change score in Dual- SeqDB (adjusted p-value  0.05). The average log2 fold change is Av S  PW RS ¼ ð2Þ assigned to the query protein and the standard deviation is stored. max ðAv S  PW Þ Queries with no remaining hits are labeled as ‘‘NA”. The standard 0–1 normalization is performed at the end, assigning a value of 1 2.4. Model validation to the highest reported fold change score, and a value of 0 to the lowest. To validate the performance of our algorithm, each individual 3. Sequence alignment against PHI-base: We downloaded PHI- PPI from the positive dataset was taken out of the predictive mod- base (accessed 2 September 2021), a dataset containing informa- els and used as query input (leave-one-out cross-validation). All tion on the role of pathogenic genes in bacteria. PHI-base assigns the PPIs recovered in the combined interactome generated were to each entry a ‘‘mutant phenotype”, depending on how the gene considered True Positives (TP), while the rest were considered deletion or mutation affects the pathogenicity of the organism. In False Negatives (FN). Subsequently, we passed each of the non- some cases, the same gene can have more than one entry, since interacting pairs of proteins as input to our predictive algorithm. it may have been measured in different experiments. We filtered In this case, any non-interacting pair recovered in the combined out those entries referring to pathogens that do not belong to the interactome was treated as a False Positive (FP), whereas the bacterial kingdom. Then, we only kept entries with mutant pheno- remaining ones were treated as True Positives (TP). Evaluation type tags that matched ‘‘unaffected pathogenicity”, ‘‘loss of metrics were calculated as: pathogenicity”, ‘‘reduced virulence”, ‘‘lethal” or ‘‘increased viru- lence (hypervirulence)”, and transformed them into numerical val- TP Precision ¼ ð3aÞ ues, 0, 0.5 or 1 as follows: ‘‘lethal” = 1, ‘‘loss of pathogenicity” = 1, TP þ FP ‘‘reduced virulence” = 0.5, ‘‘increased virulence (hypervirulence)” = TP 0.5, ‘‘unaffected pathogenicity” = 0. Recall ¼ ð3bÞ In the event that discrepant phenotypes were reported for the TP þ FN same database entry, the most abundant tag is assigned. HPIPred TN then performs a protein sequence alignment between each query Specificity ¼ ð3cÞ protein and PHI-base. Hits with a percentage of identity  40 % TN þ FP and an E-value < 10 are retained. Surviving queries are assigned FP an average PHI-base score and the mean standard deviation is Sensitiv ity ¼ ð3dÞ FP þ TN stored. ‘‘NA” labels are assigned to queries with no hits. 4. Betweenness centrality of host proteins: As suggested by the TP þ TN centrality-lethality rule [23,24], proteins that are central in the Accuracy ¼ ð3eÞ TP þ FP þ TN þ FN interactome are more likely to be essential for the organism. In this sense, betweenness centrality (BC) is a relevant centrality measure, TP  TN  FP  FN as nodes with high betweenness are located on key communica- MCC ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð3fÞ ðTP þ FPÞ  ðTP þ FNÞ  ðTN þ FPÞ  ðTN þ FNÞ tion routes and control network integrity [25]. Hence, we used BC as a proxy for protein relevance in the host. To measure protein essentiality in the Homo sapiens proteome, we calculated the BC 2.5. Software implementation score for all proteins. The Homo sapiens interactome was down- loaded from the STRING database [26] (accessed 2 September HPIPred was designed in R and wrapped up in a command-line 2021. We filtered out all PPIs with a confidence score lower than tool that allows the user to choose the host and pathogen organ- 0.9. We then used the R-package igraph [27] to build an undirected isms. The software can be used in Linux and MacOS operating sys- network, calculated the node BC score for all nodes in the graph, tems and is available at: https://github.com/SysBioUAB/hpi_ each representing a human protein, and performed a standard 0– predictor. Users can select a Uniprot proteome ID or upload the 1 normalization, being 0 the protein with the lowest BC score proteomes as custom files. They can also choose among more than and 1 with the protein with the highest score (Fig. 3B). 400 different physicochemical descriptors, select the preferred 5. Calculation of the ranked score: For each predicted PPI in the false positive rate and decide the number of models to build the combined interactome we compiled all the normalized scores consensus interactome. We also provide pre-calculated protein obtained in the previous steps (BacFITBase, DualSeqDB, PHI-base similarity values for five model organisms that can be used as scores for the pathogen proteins, and betweenness centrality and hosts, namely Homo sapiens, Mus musculus, Dario rerio, Caenorhab- DualSeqDB scores for the host proteins) and calculated an average ditis elegans, Drosophila melanogaster, for the five default physico- score with a value ranging from 0 to 1 using the following chemical descriptors. This option speeds up calculations when equation: the model organisms are used. The calculation of an interactome 6538 J. Macho Rendón, R. Rebollido-Ríos and M. Torrent Burgas Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 can take one day to compute in an average computer, while the pre-calculated proteomes can decrease the time to a few hours. 3. Results 3.1. Model prediction scores One of the main goals in PPI prediction is keeping the number of false positives low. Hence, for algorithms built to predict inter- species interactomes, it is essential to control the false positive rate (FPR). Otherwise, the number of wrongly predicted interactions would bias the interpretation of the data and undermine the efforts to validate the results in the laboratory. When predicting bacteria- human interactomes, the number of potential interactions to test is around 80 million, assuming a bacteria proteome size of 4.000 pro- teins (20.000 human  4.000 bacteria proteins). In this scenario, even with a small FPR of 0.1 %, the number of wrong predictions would be in the range of 80.000 interactions. In our method, we used synthetic proteomes, to keep these numbers low. These non-biological proteomes allowed us to make a coarse estimate Fig. 4. Model evaluation of different PPI prediction methods. Sequence-based of the FPR in large-scale predictions and control for low FPRs. methods are based on numerical descriptors, annotation-based methods are based Specifically, we set the maximum FPR to 0.0001 %, giving a theoret- on domain and GO annotations, and HPIPred is the method described here. Control results were obtained based on equivalent sequence-based methods but using ical estimate of 80 false positive interactions in a standard random numerical vectors instead of meaningful descriptors. Measurements for interactome. control, sequence- and annotation-based methods, were obtained from B. Dunham To validate our predictions, we used the leave-one-out strategy. et al.[16]. Hence, all entries in the positive dataset were used as templates to evaluate every interaction, except the one being tested. Using default settings, our method virtually achieves a precision of predicted interactions. The prediction of interactomes of relatively 100 % when combining all five models (Supplementary Informa- small sizes but with high confidence can be seen as an advantage tion, Table S1). However, this comes at the cost of a low recall, for downstream analysis, including network visualization or gene i.e., only 2 % of the true interactions are recovered, as observed ontology enrichment, but also for further experimental validation. in the ROC and PR curves (Supplementary Information, Fig. S1). We chose to further explore the predicted interactome gener- Although we achieved a low recall, the predictive power of our ated with the most restrictive filter (FPR = 0.00001 %) due to its method is similar when compared with other protein–protein pre- suitable proteome size (763 PPIs). After calculating the node diction algorithms [16], but with the advantage of controlling the betweenness centrality of the proteins in the Homo sapiens pro- FPR to very low levels (Fig. 4). In most sequence-based prediction teome and performing sequence alignment against DualSeqDB, algorithms, precision can drop to 10 % at 3 % recall and 20 % in BacFITBase, and PHI-base, we generated a ranked score for each annotation-based methods. All prediction methods, including our PPI. This, in turn, allowed us to prioritize the interactions not only own, suffer mainly from using incomplete datasets of PPIs: since according to protein similarity based on physicochemical proper- only a small fraction of the search space of PPIs has been experi- ties but also on network topology and biological properties related mentally validated, the models do not perform well for generaliza- to the infection process. We visualized the predicted interactome tion, because the penalty of removing a known PPI from the in Cytoscape [28] (Fig. 4) by representing it as a bipartite graph, positive dataset is very costly. Furthermore, for host-pathogen where host proteins (host nodes) are only connected to pathogen interactomes, only a few species of bacteria have been studied, proteins (pathogen nodes) and vice versa. The ranked score of each which also limits the available information on which these predic- predicted PPI was represented by the thickness and the color tors are based. intensity of the edges, that is, the closer a ranked score to 1, the thicker its edge and the more intense its color, meaning that this 3.2. Prediction of the complete host-pathogen interactome Homo PPI scored well on average on the biological databases. sapiens – Pseudomonas aeruginosa PAO1 In addition, we colored the nodes to represent changes in expression, derived from the BLAST search against DualSeqDB, in To test how our model performs in a real case scenario, we pre- a range from blue (downregulated) to red (upregulated). We high- dicted the Homo sapiens – Pseudomonas aeruginosa PAO1 interac- lighted some of the top scoring PPIs in the Cytoscape network to tome. The H. sapiens and P. aeruginosa proteomes contain 19.192 show how highly ranked PPIs (Fig. 5), which represent protein and 1.314 proteins, respectively, after length filtering, resulting pairs with inferred biological relevance, also appear to be impor- in 25 million putative interactions. We calculated the CCCs for tant for network integrity and connectivity. all these interactions and generated the interactomes at 5 different We then filtered out all PPIs with a ranked score lower or equal FPRs (0.1 %, 0.01 %, 0.001 %, 0.0001 %, and 0.00001 %) with five to 0.6 and divided the proteins involved in these PPIs into a patho- physicochemical descriptors (GEIM800101, GEIM800105, gen fraction (16 proteins) and a host fraction (128 proteins) of ZIMJ680101, ZIMJ680104, COSI940101), as well as their combined unique proteins, which were then used to perform gene ontology model (Table 1). As noted before, in the case of single models, the enrichment analysis with DAVID [29] (Fig. 6). The results showed allowed FPR can restrict the number of interactions from a hun- that the host fraction of proteins was enriched in biological pro- dred thousand PPIs when using the most permissive filter (0.1 %) cesses related to immune and inflammatory response, such as reg- to only several hundred PPIs with the most restrictive filter ulation of NF-KappaB activity, cellular iron homeostasis, and actin (0.00001 %). Also, the combination of different models can decrease filament bundle assembly. As for the pathogen fraction, we the number of predicted PPIs but increase the confidence in the observed an enrichment in biological terms related to amino acid 6539 J. Macho Rendón, R. Rebollido-Ríos and M. Torrent Burgas Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 Table 1 PPI sizes of the predicted interactomes by individual and combined models, at different FPRs. Predicted interactions FPR Model ZIM101 ZIM104 GEIM101 GEIM105 COSI940101 Combined 0.1 % 142.089 386.529 551.043 143.493 225.456 2.594 0.01 % 8.469 15.155 22.907 7.082 15.119 1.661 0.001 % 1.853 2.322 3.491 1.485 3.717 1.104 0.0001 % 1.035 1.085 1.439 968 1.779 894 0.00001 % 814 915 915 763 1.146 763 Fig. 5. Network representation of the Homo sapiens – Pseudomonas aeruginosa PAO1 interactome predicted by the combined models. Nodes represent proteins and edges represent predicted interactions between proteins. Host and pathogen proteins are represented by circles and squares, respectively. Nodes are colored according to the normalized expression changes (computationally derived from DualSeqDB) or black in case of missing information. Node sizes are proportional to normalized betweenness centrality and fitness for the host and pathogen proteins, respectively. Edge size and width correspond to the PPI ranked score (0-1 scale). Some of the PPIs with the highest final scores have been highlighted with colors to show how the highest-ranked PPIs from our predictive algorithm allow reconstructing of a highly connected subnetwork. The network was designed with Cytoscape. Fig. 6. Gene Ontology enrichment analysis of proteins involved in the highest scoring PPIs. Biological processes (BP) are depicted for (A) host and (B) pathogen proteins. Log fold changes for individual proteins are assigned from sequence similarity to log2 fold changes from DualSeqDB. Gene ontology analysis was performed using David and chord diagrams were drawn using the circlize package in R. 6540 J. Macho Rendón, R. Rebollido-Ríos and M. Torrent Burgas Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 and nucleotide biosynthesis, as well as folate biosynthesis, all of cape (Fig. 7) and found that, in general, the shared PPIs maintained them required for bacterial proliferation. These biosynthetic routes a certain degree of network connectivity and the proteins involved have been used as molecular targets for the development of presented a high score in terms of betweenness centrality and fit- antimicrobials, e.g., the antibiotic trimethoprim is a dihydrofolate ness for the host and the pathogen, respectively. These results sug- reductase inhibitor. gest that controlling for false positives is essential for useful predictions. Otherwise, a lot of time and resources would be 3.3. Benchmarking with interolog prediction servers wasted on experimental validation. We further compared HPIPred with two publicly available pre- dictive software called BIPS (Biana Interolog Prediction Server) [9] 4. Discussion and PredHPI [10]. Both algorithms predict putative PPIs based on interolog information. We used Homo sapiens and Pseudomonas Infectious diseases are a growing health concern worldwide, aeruginosa PAO1 proteomes as query inputs to compare the results specially due to the increase in multi-drug-resistant bacteria to our consensus interactome (Supplementary Information, [30]. These pathogenic bacteria cause prolonged hospitalizations, Table S2). In the case of PredHPI, we recovered 37.815 interactions higher costs for medical treatment, and increased mortality rates. using default parameters (30 % similarity). Surprisingly, some pro- In this sense, protein interactions between pathogenic bacteria teins had many interactions, such as immunoglobulin IGHV4-31, and their natural hosts play a key role in the infection mechanism efflux pump MacB and plasmin, with 8.730, 1.972, and 1.472 inter- and a thorough understanding of their complex interplay [31] is actions, respectively. In fact, 10 proteins account for 50 % of the required for the development of new antibiotics. Numerous exper- total number of interactions. These numbers are disproportionate imental techniques, such as yeast two-hybrid, pulldown assays, or and may represent a bias towards overrepresented entries in the co-immunoprecipitation, are currently used for the detection of positive databases used to train the predictor. Even increasing these interactions, which are collected in databases through litera- the similarity score to 40 %, we recovered 13.941 interactions, with ture mining or manual curation. However, experimental tech- 3.388, 177, and 156 interactions for IGHV4-31, MacB, and plasmin, niques are time-consuming, costly, and suffer from low respectively. Again, a single protein represents  25 % of the total specificity [32], making it unfeasible to evaluate all possible pro- number of interactions, suggesting a high bias in the predicted tein–protein interactions. interactome. Conversely, in BIPS, we obtained very few predictions Recently, computational approaches such as machine learning, using default parameters (<20), probably because the sequence homology-based methods, or structure-based methods, have similarity threshold was very restrictive (80 %). Hence, we relaxed allowed the prediction of putative PPIs, complementing experi- the filtering criteria involving identity similarity to 40 %. BIPS pre- mental techniques [33]. The main limitation of these prediction dictive tool generated an interactome consisting of 963 PPIs, a methods lies in the high rate of false positives, mainly due to the more manageable set with a size similar to our predicted interac- lack of robust databases of experimentally validated PPIs [16]. This tome. The results from BIPS and our algorithm were then com- data shortage causes the prediction methods to perform poorly in pared by creating an intersection of the predicted interactomes, terms of generalization. HPIPred has been developed as a tool that which revealed that both methods shared a total of 262 common could help to reduce the false positive rate compared to other PPIs. We represented these common PPIs as a network in Cytos- methods. To this end, HPIPred predicts putative PPIs through the Fig. 7. Network representation of the common PPIs by BIPS and HPIPred. Host and pathogen proteins are represented by circles and squares, respectively. Nodes are colored according to the normalized expression changes (computationally derived from DualSeqDB) or black in case of missing information. Node sizes are proportional to normalized betweenness centrality and fitness for the host and pathogen proteins, respectively. Edge size and width correspond to the PPI ranked score (0-1 scale). The network was designed with Cytoscape. 6541 J. Macho Rendón, R. Rebollido-Ríos and M. Torrent Burgas Computational and Structural Biotechnology Journal 20 (2022) 6534–6542 numerical encoding of proteins based on physicochemical proper- [4] Cossar PJ, Lewis PJ, McCluskey A. Protein-protein interactions as antibiotic targets: a medicinal chemistry perspective. Med Res Rev 2020;40:469–94. ties of the amino acids and integrates these predictions with bio- [5] Karimova G, Ladant D, Ullmann A. Two-hybrid systems and their usage in logically relevant data. Such data include information on the infection biology. Int J Med Microbiol 2002;292:17–25. in vivo relevance of bacterial genes to the infection process, gene [6] Gagarinova A, Phanse S, Cygler M, Babu M. Insights from protein-protein expression changes in vitro and in vivo, as well as topology infor- interaction studies on bacterial pathogenesis. Expert Rev Proteomics 2017;14:779–97. mation that allows highlighting the importance of central hubs. [7] Jean Beltran PM, Federspiel JD, Sheng X, Cristea IM. Proteomics and integrative By using the Homo sapiens and the Pseudomonas aeruginosa PAO1 omic approaches for understanding host-pathogen interactions and infectious proteomes as input to our prediction tool, we generated 763 diseases. Mol Syst Biol 2017;13:922. [8] Gómez Borrego J, Torrent Burgas M. Analysis of Host-Bacteria Protein host-pathogen interactions displaying a highly connected network. Interactions Reveals Conserved Domains and Motifs That Mediate We expect that our prediction tool will provide a more realistic Fundamental Infection Pathways. International Journal of Molecular Sciences picture of host-pathogen interactomes and help pave the way for 2022;23(19):11489. [9] Garcia-Garcia J, Schleker S, Klein-Seetharaman J, Oliva B. BIPS: BIANA Interolog the prioritization of PPIs that can be explored as potential targets Prediction Server. A tool for protein-protein interaction inference. Nucl Acids for the development of new antibacterial drugs. Res 2012;40:W147–51. [10] Loaiza CD, Kaundal R. PredHPI: an integrated web server platform for the detection and visualization of host-pathogen interactions using sequence- CRediT authorship contribution statement based methods. Bioinformatics 2021;37:622–4. [11] Das D, Krishnan SR, Bulusu G, Roy A. in 2019 IEEE International Conference on Javier Macho Rendón: Data curation, Formal analysis, Method- Bioinformatics and Biomedicine (BIBM) 935-938 (2019). [12] Zhang SB, Tang QR. Protein-protein interaction inference based on semantic ology, Software, Writing – original draft. Rocio Rebollido-Ríos: similarity of Gene Ontology terms. J Theor Biol 2016;401:30–7. Data curation, Formal analysis, Methodology, Writing – review & [13] Mariano R, Wuchty S. Structure-based prediction of host-pathogen protein editing. Marc Torrent Burgas: Supervision, Writing – review & interactions. Curr Opin Struct Biol 2017;44:119–24. [14] Kaundal R, Loaiza CD, Duhan N, Flann N. deepHPI: a comprehensive deep editing, Funding acquisition. learning platform for accurate prediction and visualization of host-pathogen protein-protein interactions. Brief Bioinform 2022;23. Declaration of Competing Interest [15] Li F, Zhu F, Ling X, Liu Q. protein interaction network reconstruction through ensemble deep learning with attention mechanism. Front Bioeng Biotechnol 2020;8:390. The authors declare that they have no known competing finan- [16] Dunham B, Ganapathiraju MK. Benchmark evaluation of protein-protein cial interests or personal relationships that could have appeared interaction prediction algorithms. Molecules 2021;27. [17] Durmus Tekir, S. et al. PHISTO: pathogen-host interaction search tool. to influence the work reported in this paper. Bioinformatics 29, 1357-1358 (2013). [18] Tiessen A, Perez-Rodriguez P, Delaye-Arredondo LJ. Mathematical modeling Acknowledgments and comparison of protein size distribution in different plant, animal, fungal and microbial species reveals a negative correlation between protein size and protein number, thus providing insight into the evolution of proteomes. BMC This study was funded by a Research Grant 2022 of the European Res Notes 2012;5:85. Society of Clinical Microbiology and Infectious Diseases (ESCMID) [19] Kawashima S, Kanehisa M. AAindex: amino acid index database. Nucl Acids and the Spanish Ministerio de Ciencia e Innovación (SAF2017- Res 2000;28:374. [20] Rendon JM, Lang B, Tartaglia GG, Burgas MT. BacFITBase: a database to assess 82158-R, PDC2021-121544-I00 funded by MCIN/AEI the relevance of bacterial genes during host infection. Nucl Acids Res 2020;48: /10.13039/501100011033 and European Union Next Genera- D511–6. tionEU/ PRTR, and project PID2020-114627RB-I00 funded by [21] Macho Rendon J, Lang B, Ramos Llorens M, Gaetano Tartaglia G, Torrent Burgas M. DualSeqDB: the host-pathogen dual RNA sequencing database for infection MCIN/AEI /10.13039/501100011033), all to MT. R.R.-R was a recip- processes. Nucl Acids Res 2021;49:D687–93. ient of an INCOMED Marie Curie Fellowship at the time. [22] Urban M et al. PHI-base: the pathogen-host interactions database. Nucl Acids Res 2020;48:D613–20. [23] Jeong H, Mason SP, Barabasi AL, Oltvai ZN. Lethality and centrality in protein Author contributions networks. Nature 2001;411:41–2. [24] Crua Asensio N, Munoz Giner E, de Groot NS, Torrent Burgas M. Centrality in MT designed, directed, obtained funding for, and coordinated the host-pathogen interactome is associated with pathogen fitness during infection. Nat Commun 2017;8:14092. the study. JM and RR-R performed the experiments. JM wrote an [25] Ashtiani M et al. A systematic survey of centrality measures for protein- initial version of the paper, subsequently edited by MT. All authors protein interaction networks. BMC Syst Biol 2018;12:80. contributed to editing and revising the final version of the [26] Szklarczyk D et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/ manuscript. measurement sets. Nucl Acids Res 2021;49:D605–12. [27] Csardi G, Nepusz T. The igraph software package for complex network research. InterJournal. Complex Systems 2006;1695. Appendix A. Supplementary data [28] Shannon P et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res 2003;13:2498–504. Supplementary data to this article can be found online at [29] Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc https://doi.org/10.1016/j.csbj.2022.11.026. 2009;4:44–57. [30] Antimicrobial Resistance C. Global burden of bacterial antimicrobial resistance References in 2019: a systematic analysis. Lancet 2022;399:629–55. [31] de Groot NS, Torrent Burgas M. A coordinated response at the transcriptome [1] Ribet D, Cossart P. How bacterial pathogens colonize their hosts and invade and interactome level is required to ensure uropathogenic escherichia coli deeper tissues. Microbes Infect 2015;17:173–83. survival during bacteremia. Microorganisms 2019;7. [2] Crua Asensio N, Macho Rendon J, Torrent Burgas M. Time-resolved [32] Zhou M, Li Q, Wang R. Current experimental methods for characterizing transcriptional profiling of epithelial cells infected by intracellular protein-protein interactions. ChemMedChem 2016;11:738–56. Acinetobacter baumannii. Microorganisms 2021;9. [33] Keskin O, Tuncbag N, Gursoy A. Predicting protein-protein interactions from [3] de Groot NS, Torrent Burgas M. Bacteria use structural imperfect mimicry to the molecular to the proteome level. Chem Rev 2016;116:4884–909. hijack the host interactome. PLoS Comput Biol 2020;16:e1008395. 6542