Abstract

Motivation: Analysis of protein–protein interaction networks (PPINs) at the system level has become increasingly important in understanding biological processes. Comparison of the interactomes of different species not only provides a better understanding of species evolution but also helps with detecting conserved functional components and in function prediction.

Method and Results: Here we report a PPIN alignment method, called PINALOG, which combines information from protein sequence, function and network topology. Alignment of human and yeast PPINs reveals several conserved subnetworks between them that participate in similar biological processes, notably the proteasome and transcription related processes. PINALOG has been tested for its power in protein complex prediction as well as function prediction. Comparison with PSI-BLAST in predicting protein function in the twilight zone also shows that PINALOG is valuable in predicting protein function.

Availability and implementation: The PINALOG web-server is freely available fromhttp://www.sbg.bio.ic.ac.uk/~pinalog. The PINALOG program and associated data are available from the Download section of the web-server.

Contact:m.sternberg@imperial.ac.uk

Supplementary information:Supplementary data are available atBioinformatics online.

1 INTRODUCTION

In recent years, there has been a substantial increase in the quantity of protein–protein interactions (PPIs) being detected by experimental methods. The analysis of PPI networks (PPINs) at the system level has become increasingly important in understanding biological processes. Comparison of the interactomes of different species not only provides a better understanding of species evolution but also helps with detecting conserved functional components and in function prediction.

Several methods have been developed to align interactomes, both globally and locally, most of which use sequence similarity or network topology or both in establishing the equivalence. PathBLAST (Kelleyet al., 2003) marked the first local alignment method and was similar in principle to the BLAST search algorithm. Other local alignment methods include: NetworkBLAST (Sharanet al., 2005), which is the next generation of PathBLAST; MaWISH (Koyutürket al., 2006), which adopts the evolutionary models of match, mismatch and deletion of the proteins; Graemlin 1.0 (Flannicket al., 2006), which infers an alignment from network modules; the Bayesian method (Berg and Lassig, 2006); the match and split algorithm (Narayanan and Karp, 2007); and Phunkee (Cooteset al., 2007), which aligns proteins based on sequence and context in the networks. The global alignment of networks proves to be more challenging due to the complexity and scale of the problem. Graemlin 2.0 (Flannicket al., 2009) formulates a model for protein duplication, deletion and mutation and aligns the network progressively using a hill-climbing algorithm. The Markov-random field-based method (Bandyopadhyay, 2006) and IsoRank (Singhet al., 2008) solve the problem by eigenvalue-based methods. GRAAL (Kuchaievet al., 2010) uses only the network topological similarity to align the networks using a seed-and-extend method for graph alignment.

In many applications of PPIN alignment, the underlying objective is to equivalence proteins with related function and interacting partners. However, in most alignment methods only the sequence data or network topology are used to align the input networks despite the availability of other sources of information such as function annotation or gene expression. While sequences are informative in elucidating the orthologous relationships of proteins across species, they do not necessarily indicate functional similarity. The omission of function information could therefore result in many pairs of equivalenced proteins having little or no similarity of function. This in turn makes it less accurate in detecting conserved functional modules or predicting protein function. Recently,Ali and Deane (2009) have introduced functional similarity of proteins into a local alignment method. They used functional similarity in combination with sequence similarity in the local alignment match and split algorithm (Narayanan and Karp, 2007) and identified similar functional subnetworks. The recently developed method MI-GRAAL (Kuchaiev and Przulj, 2011) presents a global alignment algorithm in which different information can be incorporated such as topological features, sequence similarity and functional similarity. However, there has been no assessment on how the inclusion of functional similarity might influence the resulting alignments and any subsequent applications.

Here we describe PINALOG, a global network alignment algorithm. PINALOG forms the alignment between two PPINs based on the similarities of protein sequence and the protein function between the two networks. Functional similarity is formalized using GO (gene ontology) annotations. The main contributions in PINALOG are the use of communities in the networks to identify seed protein pairs and the scoring schemes used in the extension steps to include the neighbourhood similarity of mapped protein pairs. The benchmarking of PINALOG in comparison with other alignment methods shows that PINALOG obtains a good balance in including these features as it generates an alignment with consistently high number of conserved interactions, homologous aligned pairs and functionally similar protein pairs. Direct applications of PINALOG include protein complex prediction and protein function prediction. Complex prediction is performed by direct inheritance of protein complex data from a known species to an unknown one. PINALOG has been tested for complex prediction power using a human-yeast alignment. Protein function prediction is via direct transfer to an un-annotated protein from an annotated protein. The power of our function prediction is assessed in a cross-validated study.

2 METHODS

The design of the PINALOG alignment method is based on the observation that proteins are not uniformly distributed throughout PPINs. Instead there are some proteins that form well connected subnetworks. Furthermore, it has been shown byBrohee and van Helden (2006) that protein complexes and functional modules form highly connected components in the PPINs. Therefore, it would be more reliable and efficient to align two PPINs by first finding highly similar protein pairs (seed protein pairs) from these highly connected protein subnetworks (referred to as communities from now on) in the networks and then extending the alignment to other proteins in the neighbourhoods of seed protein pairs.

LetA andB be the PPINs for species A and B.a_i andb_j are thei-th andj-th proteins in PPIN A and B, respectively.C_i^A is thei-th community in PPINA, andC_j^B is thej-th community inB. PINALOG alignsA andB in three steps as summarized inFigure 1: community detection, community mapping and extension mapping.

Step 1: Community detection of input networks using CFinder

Communities in biological networks such as PPINs often indicate functional groupings of proteins in the network (Adamcseket al., 2006;Lewiset al., 2010;Song and Singh, 2009). There are several methods to detect communities within a network such as the minimum-cut method, hierarchical clustering or the Girvan–Newman algorithm. In PINALOG, we employ CFinder (Pallaet al., 2005), a clique percolation method, which is capable of detecting overlapping communities. Denote ak-clique as a subgraph of the network composed ofk proteins where all pairs of proteins interact. Then ak-clique community is defined as union of allk-cliques that are reachable from each other through adjacentk-cliques (two adjacentk-cliques havek−1 proteins in common). CFinder constructs communities by merging adjacentk-cliques (seeSupplementary Material for more details).

Node scoring scheme

The similarity between two proteins is a combination of sequence and functional similarity in the seed identification and with additional topological similarity in the extension mapping process. The sequence similarity of two proteinsa_i andb_j is calculated based on their BLAST bit score as

(1)

S(a_i,b_j) is the BLAST bit score value when aligninga_i andb_j. Here, only pairs of proteins with anE-value <10⁻⁵ are used to calculate sequence similarity.

The functional similarity of two proteins annotated with GO terms is calculated by the method defined bySchlickeret al. (2006). Given two proteinsa_i andb_j with two sets of associated GO termsg(a_i)={t₁,t₂,…,t_k} andg(b_j)={t′₁,t′₂,…,t′_l}, the protein functional similarity is calculated based on the semantic similarity (Lordet al., 2003) of terms ing(a_i) andg(b_j). Semantic similarity depends on the rarity of the terms in the annotation database as well as their distance in the ontology. The semantic similarity between each term ing(a_i) and each term ing(b_j) is calculated, resulting in ak×l table of semantic similarities. Let rowScore and colScore be the average of the row maxima and the average of column maxima in the table, respectively (seeSupplementary Fig. SI-1), the Schlicker's similarity betweena_i andb_j is then calculated as

(2)

The similarity of two proteins in Step (2) is defined as a linear combination of functional similarity and sequence similarity

(3)

In PINALOG, θ provides a relative weighting between sequence and functional similarity. The value of θ is automatically decided based on the number of reciprocal best BLAST hits in the protein sequences between two input species, reflecting the closeness of two species and thus the contribution of sequence similarity in the overall score. Details of how θ is calculated are given inSupplementary Material.

Step 2: Community mapping to obtain seed protein pairs

Community mapping is determined in two steps. The first step is to determine the highest score when equivalencing proteins in each community in species A with each community in species B. To this end, we define the similarity between two communities as the scoreF(C_i^A,C_j^B), the sum of similarities between protein pairs in the optimal equivalence (OptMap) of proteins in communityC_i^A in species A with proteins inC_j^B in species B using the Hungarian method (Kuhn, 2005). The optimal equivalence is the mapping where the sum of protein pair similarities is largest. The Hungarian algorithm is a combinatorial optimization algorithm solving the assignment problem in polynomial time (seeSupplementary Material).

(4)

The first step produces a community similarity matrix. This matrix is used in the second step to obtain the best equivalence of communities by maximizing the community scoring function:

(5)

The maximization ofF(core) is achieved by Hungarian method based on the community similarity matrix. For each pair of mapped communities, the Hungarian mapping of proteins between them are obtained and added to the list of seed protein pairs. A filtering step is used to retain only the top 15% pairs of mapped proteins as the seeds for extension, known as the core equivalences.

Step 3: Extension mapping

Extension mapping is performed by considering neighbours of proteins in the core equivalence as candidates for adding to the alignment. In addition to protein sequence and functional similarity, topological similarity in the PPINs is included in the form of neighbourhood similarity. LetN(a_i) andN(b_j) be the set of all first neighbours (proteins separated by one interaction) and second neighbours (proteins separated by two interactions) ofa_i inA andb_j in B. Letd(a_k,a_l) denote the distance betweena_k anda_l in a network. The similarity betweena_i andb_j in extension mapping is then defined as

(6)

By awarding the candidate pairs with a proportion of the score of neighbouring aligned pairs, PINALOG aims at mapping more neighbouring protein pairs. The Hungarian method is used to find the optimal equivalence of candidates. These candidates are then added to the equivalences in the core. The process is repeated until no more pairs are added (seeSupplementary Fig. SI-1). As CFinder identifies communities that can be overlap, community mapping can result in many-to-many mappings in seed protein pairs, thus in the final alignment.

3 RESULTS

4 PROTEIN COMPLEX PREDICTION

One of the applications of network alignment is to use the available information on the protein complex in one species to predict the protein complex components in another species. To benchmark the protein complex prediction of PINALOG, we performed a comparison of PINALOG with network-based complex prediction methods, including MCODE (Bader and Hogue, 2003), MCL (Dongen, 2000), network alignment methods including IsoRank, MI-GRAAL and BLAST. Comparisons were made on the gold-standard set of yeast complexes MIPS CYDG (Güldeneret al., 2005) comprising of 214 protein complexes. Predictions of the yeast complexes were made by PINALOG, IsoRank, MI-GRAAL and BLAST by transferring the protein complex information from human to yeast, with the human protein complexes information obtained from the MIPS CORUM (Rueppet al., 2008) database.

An assessment of prediction methods was performed using a series of functions based on the overlapping score (OS, seeSupplementary Material). The predictionp is said to match with complexm if OS(p,m) is >0.2. Given the scores of matching between predictions and known complexes, measures to assess the quality of prediction include: (i) PrecisionP is the fraction of predicted complexesp matching known complexesm; (ii) RecallR is the fraction of known complexesm matching predicted complexesp; and (iii)F-measureF₁ is the harmonic mean of precision and recall, providing the overall performance of the prediction methods.

Figure 2 presents the precision, recall andF-measure of yeast protein complex prediction by PINALOG-A. PINALOG outperforms MCODE, MCL, IsoRank and MI-GRAAL in all three measures, with 37% recall at 21% precision and anF-measure of 26%. The performance of PINALOG is better than BLAST in precision but worse in recall. However, in theF-measure, which summarizes the overall performance, these two methods perform equally well.

It should be noted that the precision and recall are not high because the benchmark protein complexes are those identified by experiments that do not cover all protein complexes in yeast. Thus a high number of false positives (FPs) are expected. These FPs are candidates for protein complexes that could be studied experimentally. To assess whether these FPs are likely to be true protein complexes or not, we have examined the predicted clusters that do not match with known complexes. Using GOTermFinder (Boyleet al., 2004), we have identified significant GO term enrichment in many unmatched predicted clusters. For example, the 9 proteins of cluster 13 are enriched with terms related tochromatin modification,chromatin organization, andchromatin assembly and disassembly. Cluster 2, comprising 23 proteins, is enriched with terms relating tooxidation–reduction process. Another example is cluster 12 with 8 proteins that are related toM phase duringcell division. The GOTermFinder analysis indicates that the putative FPs in the PINALOG prediction may well be true protein complexes.

BLAST performed as well as PINALOG and better than other single network-based methods. However, when inspecting the overlap in the predictions made by PINALOG and BLAST, we found that the overlap of the two predicted sets of clusters is not substantial, with only 50% of the predicted clusters matched. This suggests that PINALOG and BLAST could be used complementarily to predict protein complexes. This would provide a better coverage of the possible protein complexes. BLAST predicts clusters based on protein sequence, thus the prediction space is limited to sequence related proteins. In contrast, PINALOG predictions are made based on a combination of sequence, function and network topology, therefore the prediction space is extended to functionally related protein clusters that are reinforced by conserved interactions. We have combined the predictions made by PINALOG and BLAST into one predictor PINALOG+BLAST. This combination not only boosted the precision but also the recall rate, yielding a higherF-measure as compared with the individual methods (Fig. 2). We suggest that PINALOG and BLAST are used combination to help predicting protein complex more accurately.

5 BENCHMARKING OF PROTEIN FUNCTION PREDICTION

Another application of PINALOG is to predict the function of proteins by inheriting the annotation available of the aligned protein from the other species. PINALOG aims at providing function prediction (i) when there is no sequence homology; (ii) when there is a sequence homology between the un-annotated protein and the aligned annotated protein in the other species but the per cent sequence identity is low (typically <30%) and thus direct functional transfer can lead to misleading annotations (Tian and Skolnick, 2003;Toddet al., 2001). PINALOG is compared with the widely-used sequence-based method PSI-BLAST and to IsoRank.

A dataset of proteins to test the accuracy of function prediction was established with 415 GO annotated human proteins. These proteins had low per cent identity (<30%) to the closest annotated PSI-BLAST hit in the UniprotKB database. A 100-fold cross validation was performed. In each run of the cross validation, a set of proteins from the test set and had their GO terms hidden. The test proteins which were aligned with an annotated protein in the other species were identified and the function transferred. Over the 100 runs, a total of 169 proteins were equivalenced and their function predicted. Functions of mapped proteins in yeast represented by GO terms in two categories BP (biological process) and MF (molecular function) were directly transferred to the human proteins in the test set.

The results were analyzed in terms of precision and recall, where precision is defined as and recall as; tp, fp and fn being the number of true positives, false positives and false negatives of the predictions made. Our method is a binary predictor without a variable cut-off parameter and thus no precision/recall curve can be produced.Table 2 shows that PINALOG outperforms PSI-BLAST prediction with a superior recall rate at similar level of precision with BP category. The McNemar test (McNemar, 1947) for the statistical significance of the difference in the performance was used based on the number of misclassifications in each method. The test indicates that PINALOG predictions are significantly different from PSI-BLAST in BP terms at the level ofP=0.001 significance level. PINALOG predicts more BP terms with a recall twice that of PSI-BLAST (14% versus 7%) at the same level of precision (~28%). For the MF category, PINALOG and PSI-BLAST share very similar levels of recall (~28%) while the precision of PSI-BLAST is slightly better (43% versus 47%). The McNemar test suggests no significant difference between them at theP=0.001 significance level. On the other hand, PINALOG outperforms IsoRank in both BP (P=0.001) and MF (P=0.001) categories. For example, in the BP category, PINALOG has almost twice the recall compared with IsoRank (14% versus 8%) at higher precision (28% versus 20%,Table 2). To summarize, in the challenging area where sequence similarity does not contribute substantially to the prediction of protein function, PINALOG enhances the ability to predict function of these proteins.

PINALOG was used to predict functions for un-annotated human proteins without any PSI-BLAST hit in the human-yeast alignment. For the 60 such human proteins, PINALOG mapped 14 of them to the proteins in yeast and made predictions for human proteins from the yeast counterparts. The function of only nine proteins was predicted (and provided inSupplementary Table SI-2 inSupplementary Material) because the remaining five are mapped to proteins in yeast that have annotation only with very general GO terms. We compared our predictions to those from ffPred (Lobleyet al., 2007), a web-server for function prediction of orphan or un-annotated proteins. Out of these nine proteins, ffPred provides functions for six proteins, three of which have very general GO terms, the rest of the predictions agreeing with our predictions. We have also run these nine proteins through the structure prediction server Phyre2 (Kelley and Sternberg, 2009) which detects templates which can be remote homologues not identifiable by PSI-BLAST. For six of the PINALOG predictions, the functional descriptions of these templates proteins match with or are similar to the PINALOG predictions. This strongly suggests that the PINALOG predictions are reliable and can help guide experimentalists to identify protein function.

6 DISCUSSION

Although many methods have been developed to align PPINs, comparing their performance remains difficult due to the absence of an unambiguous correct alignment, the complexity of the networks and the differing aims of the alignment methods. For local alignment methods, the difficulty is less severe as some local alignments can be evaluated by their agreement with known protein complexes. However, this approach is less useful to evaluate global alignment methods. Often assessment is performed based on the correctness of mapping orthologous protein pairs. However, this assessment measure would always penalize equivalences based on network topology when this conflicts with equivalences between orthologous pairs. Therefore this assessment method is not ideal in evaluating the performance of network alignment algorithms. Accordingly, we have also used a series of other measures.

PINALOG is flexible in allowing the use of sequence only alignment or sequence-function alignment depending on the requirement of the user. For the sequence-function alignment, parameters are automatically calculated from the input species and this will help non-experts. When aligning two species where one or both species are poorly annotated, alignment might be biased towards aligning well-annotated proteins whose functions are similar. Then it is advisable that sequence and network topology only are used to align the networks to avoid bias. We are developing a version of PINALOG to perform alignment of PPINs from multiple species. A web-server for PINALOG is available allowing users to upload necessary information and receive alignment results by email and by a link to the results file (http://www.sbg.bio.ic.ac.uk/~pinalog/). The computing time of the alignment process depends on the size of the input networks. The typical computing time is <24 h with the longest run in our assessment being the human-yeast alignment which takes 24 h on a computer with 2.8 GHz processor with 8 GB memory (seeSupplementary Material). Given the importance of the resulting alignment, this time is acceptable. The code for PINALOG is available from the website, Download section.

Supplementary Material

ACKNOWLEDGEMENT

We would like to thank Wellcome Trust for funding this project and Dr Mark Wass for valuable advice on function prediction.

Funding: Wellcome Trust

Conflict of Interest: MJES is a director and shareholder in Equinox Pharma Ltd which markets software and services in the bioinformatics and chemoinformatics area.

REFERENCES

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.