Currently, a great diversity of protein similarity measures are used to infer functions. They include sequence homology [8], [9], [10], phylogenetic ancestry [11], [12] substrate similarity [13], co-expression [14], [15], physical interaction [15], [16], [17], genetic interaction [18], [19] or analogies of sequence [20], [21] or structure motifs [22]–[27]. Some methods compare divergent, aligned proteins to spot discriminating residues that suggest functional signatures in sequences (EFICAz2) [20] or in structures (FLORA) [23]. Since relatively few of all sequence or structure variations are necessarily functionally relevant [28], [29], other methods focus on just a few but presumably key residues. For example, residues could be taken from concave protein regions since these are more likely to be functionally important (pevoSoar) [22]. More generally, residues can be taken from putative functional sites in their relative three-dimensional (3D) configuration to create 3D templates: composite structural motifs of a few amino acids that directly mediate function. Experimentally validated 3D templates are available but they do not cover all of functional space [30], however, and their simple geometric matches can be spurious [31], [32]. Profunc [33] mitigates these problems by using multiple template sources, including enzyme active sites gathered from the Catalytic Site Atlas [30], ligand and DNA binding sites [34], and automatically generated triads of amino acids called ‘reverse templates’ [33]. None of these methods begin with a dedicated approach to identify protein functional sites and their key determinants.

Here, we focus on 3D templates that embody evolutionary information. Our rationale is that functionally important residues may often be distinguished by sequence variations that correlate tightly with evolutionary divergence [35], [36] and form spatial clusters in structures [37], [38], [39]. Such clusters then suggest binding sites or catalytic sites on the surface, and allosteric pathways internally [35], [37], [40], [41], [42], [43] which, in turn, efficiently guide experiments to block, separate, rewire, or mimic function [44]–[54]. Building on these computational and experimental studies that demonstrate evolutionary identification of functional determinants, our approach ranks the relative evolutionary importance of every residue in a protein sequence with the Evolutionary Trace [35], [55] (ET), and then selects the six most important and clustered surface residues to define a 3D template. The geometric matches of these evolutionary templates in other protein structures at sites that are themselves evolutionarily important then define Evolutionary Trace Annotation (ETA) annotations [31], [56]. So far, ETA annotations have been shown to be functionally specific (positive predictive values above 90%) in enzymes and non-enzymes alike [57], but their functional resolution and coverage are limited. For example, enzyme predictions, are at the third rather than the fourth EC level and coverages range from 40 to 70 percent.

We hypothesize that we can improve performance by basing annotations not only on direct matches to a protein of interest, but on all ETA matches detected across the proteome. This follows studies that systematically pool together multiple functionally-relevant matches between proteins [58]. One such approach examines a protein-protein interaction neighborhood to choose predictions that maximize the observed versus expected frequency of a function in a local subnetwork [59]. An extension further factors the network's topological weights into the prediction [17]. Another approach employs a probabilistic analysis of the functional neighborhood defined by sequence alignment bit-scores [60]. More generally, a network in which each node is a protein and each edge is a pairwise match can be constructed and then analyzed with clique or module detection algorithms to increase the amount of information involved in a prediction. For example, Mcode uses a greedy algorithm to grow clusters from a seed node [61]. SpectralMod iteratively cuts edges until only dense clusters remain [62]. CFinder groups together tightly connected cliques into clusters [63]. Once a cluster is identified, enriched functions in that cluster are propagated to unannotated members of the cluster [64]. There is, however, no single best clustering method for all topologies, and clustering approaches can be less accurate than local methods [65].

Global graph theoretic tools may further improve performance by including the entire network topology into annotation transfer [64]. One such approach optimizes the labeling of proteins with unknown function in order to maximize connections between proteins with the same function [16]. Other methods use Markov random fields and assume that the function of every protein in the network is dependant, in probabilistic terms, only on its direct neighbors [66]. Lastly, flow or diffusion-based methods have been proposed that are able to take advantage of both global structure and local similarity. FunctionalFlow is an iterative algorithm that simulates the flow of liquid through a network for the purposes of functional annotation [67]. These methods can also be directly optimized and have been applied to prediction of protein-protein interactions [68] as well as prediction of Gene Ontology (GO) terms [69], [70].

Next, once the network has been constructed, a diffusion mechanism lets annotation data flow from node to node so that the poorly annotated regions draw information from the richly annotated ones.

Formally, let a graph have n nodes, each depicting a protein, such that a number p≤n carry labels y_i that represent specific functions, where the index I ranges from 1 to n. ETA matches provide all-against-all linkage data on by averaging the ET score and RMSD of the template match to define edges in a symmetric adjacency matrix, w. If nodes I and j have a mutual ETA template match (Figure 1A), then they are linked by edges so that the adjacency matrix entry is set to w_ij>0, otherwise w_ij=0 if there is no known or significant similarity. The problem is then to infer function for the remaining q=n - p nodes from the nodes with known labels y and the network's connectivity. This frames an optimization problem between two conditions: (i) that the function of neighboring nodes be similar (smoothness condition), and (ii) that the final label of a node be consistent with its initial label if its function was known (initial condition). Both conditions can be most simply represented by a quadratic “cost” function [71]:where the functional label y_i is set to either 1 if node I has a function y, to −1 if it does not have that function, and to 0 if there is no evidence either way (unlabeled node); and where −1≤f_i≤1 are the predicted functional labels of the output prediction vector. Finally, the coupling parameter α is analogous to a diffusion coefficient that balances smoothing with loss of the initial labels. This formulation is empirical, but it is convenient because solving f={f₁,…,f_n} is equivalent to solving a sparse linear system with the graph diffusion kernel, viz. , where is the Laplacian matrix, with the link weight matrix and , . In the more limited context of binary protein function classification the computational efficiency of this approach surpassed semi-definite programming (SDP) and support vector machines (SVM) [72].

Here, in order to extend this method to multi-label classification we need to account for the bias introduced by different label frequencies. To do this, we introduce a prediction z-score, defined by the solution as , where denotes the average and the standard deviation evaluated on all unlabeled nodes; z measures the positive deviation from the expected random mean in standard deviation units and thus eliminates any absolute bias due to initial conditions in y, which then allows a meaningful comparison between multiple functional classes of variable size within the network. After diffusing every function to every node, we use the function with the highest z-score as our prediction. In this way, we add a new step to the ETA process, coupling relevant and non-trivial edge detection with global function propagation.

We have benchmarked this method on two distinct test sets. First, we perform a leave-one-out test on the FLORA [23] test set in which we attempt to predict the function of a given protein using matches to the remaining proteins in the test set. Second, we have assembled a test set of Protein Structure Initiative [3] (PSI) proteins, and attempt to predict their function by matching them against a representative subset of the Protein Data Bank [73] (PDB). Finally, we examine the effect of lowering availability sequence similarity, an important consideration because the PSI often contains structures with few homologs [74].

Compared to other state-of-the-art approaches, confidence values proved better at sorting unreliable predictions, and in turn this improved annotation accuracy at the third and at the fourth, and most specific, EC levels. These results are general since they apply across all types of enzymes, and they are accurate since false positive rates decrease substantially—between 2- to 5-fold. The many predictions on unannotated proteins demonstrate the benefits of the repeated use of evolutionary information and its integration with structural information over the structural proteome.

In order to identify the various sources of information that improve annotation we compared the impact that negative information has on diffusion on identical PDB 90 networks and Structural Genomics testsets. Accuracy-coverage curves with negative labels (Figure S4A, red) have a substantial accuracy advantage over the same curves without negative labels (purple). Without the negative labels, accuracy falls: at 50% coverage it drops by 16% and 10.7% for 3 and 4 EC digit predictions, respectively. Thus, the -1 entries in the y vector, which indicate the knowledge that a protein does not perform a specific enzymatic function, contribute significantly to accuracy. By contrast, many annotation methods, for instance the nearest-neighbor and BLAST approaches we benchmark against, do not make use of this information. Hence, knowledge of which proteins in the network lack a particular function is critical for function prediction with diffusion, and may contribute to accuracy advantages over other methods.

Additionally, in order to assess the contribution of distant positive labels to annotation, we examined the shortest path lengths between proteins with the same and differing functions. Figure S4C shows a stacked histogram comparing the lengths of shortest paths in the network between nodes with a correct prediction and nodes in the dataset, segmented by the confidence z-score. Proteins with the same function (blue) tend to have shorter distances between them than proteins with different function (orange), indicating that functions generally cluster in the network. However the distributions have long tails, especially for predictions with a z score less than 3, so that in a number of instances proteins with the same function can be quite distant. Based on our accuracies, the diffusion process is presumably able to connect these distant proteins. Therefore, both negative labels and positive labels distant by 10 or more nodes are additional information sources that contribute to more accurate predictions in ETA Network Diffusion.

Strikingly, these results rely on the large-scale comparison of just six evolutionarily important template residues, chosen protein by protein. The accuracy of the network shows that these residues effectively capture the determinants of protein function. This in essence, validates on a large scale the notion that ET analysis identifies key functional residues—consistent with the conclusions of many experimental case studies [36]. Notably, as this study draws from many previous ideas, such as 3D templates [57], [76] evolutionary importance [35], functional site analysis [88], molecular determinants of function [23], [20] and network analysis [64], [69], [89], [90] it combines them uniquely by repeatedly relying on evolution at each step of the annotation process.

First, the 3D template residues are selected for their evolutionary importance measured by ET, and for their structural clustering. This local structural motif defined by evolution obviates the need for any prior knowledge, or assumptions, about the nature and determinants of function. This is an advantage since compared to the size of the proteome, there are relatively few proteins with reliable data on the molecular mechanisms underlying their function and specificity, as may be available from the catalytic site atlas [34]. Likewise, these 3D templates also replace searches for structural features, such as clefts, cavities or depressions, which are suggestive but rarely sufficient [88].

Second, the selection of which 3D template matches are functionally relevant also relies on evolution. Out of the profusion of purely geometric matches between a template and protein structures, only those that involve residues with evolutionary importance similar to the template residues are retained. Every accepted match, and therefore every edge in the network, indicates reciprocal similarities of evolutionary constraint and structural context, which raises the likelihood of a true functional similarity.

It follows, third, that diffusion over a network defined by these evolutionary template matches disseminates evolution-guided inferences over the structural proteome. The correlation between the confidence z-score associated with every diffused function and the reliability of annotations confirm that these three hierarchical types of evolutionary inferences—meaning the 3D template, the match, and the diffusion—are all well founded: evolutionary analysis thus dramatically narrows the search for the essential determinants of a protein's function and for their comparison.

The global network approach also has many intrinsic advantages. It removes the heuristic aspect of the ETA voting approach, [57] it enables global and formal integration of information over the entire structural proteome, and, as a future direction, it prepares the integration of ETA information with many other types of functionally relevant protein similarity, since the latter usually come in the form of pairwise relationships are therefore well suited for network representation [64], [69]. Specifically, the diffusion process is non-local and draws information from all of the functional labels in the network, not just those from direct matches. As a result, it extends prediction coverage compared to strictly local techniques. This is illustrated, for example, by the gene PHO147 in Pyrococcus horikoshii (PDB 2dz9 chain A), as shown in Figure 6. This protein matches solely unnannotated proteins in a well-connected cluster, so both ETA and nearest neighbors can make no prediction. Network diffusion, however, enables more distant annotations to inform the annotation of this node, leading to a correct fourth EC digit prediction.

The computation of a confidence z-score that correlates with prediction accuracy is another contribution of this work. Together with the global diffusion process, it enables unbiased consideration of all possible functions, establishes an objective criterion for selecting the best candidate, and attaches a confidence value to it. As a result, predictions can be stratified by the z-score, yielding the accuracy versus coverage receiver-operator curves, shown in Figures 3C and 3D, that remain close to 100% accuracy over a large coverage. Predictions with a z-score above 2 are over 99% accurate and over a z-score of 1 they are 98% accurate. Conversely, on the FLORA set, the vast majority of false positives also had the lowest confidences (z<−0.05) (Figure 2). The z-score is therefore an adequate marker of confidence with which to recognize unreliable predictions that otherwise would become false positives. Overall, we see from 2 to 5 fold reductions in false positives when compared to ETA, FLORA, nearest neighbors or BLAST.

Collectively these improvements yield confident predictions at the 4^th EC level, which identifies precise substrates in many cases. For example, the predicted EC annotation for gene PHO147 in Pyrococcus horikoshii (PDB structure 2dz9A) is biotin—acetyl-CoA-carboxylase ligase (EC 6.3.4.15). This function indicates the substrates ATP, biotin and acetyl-CoA-carboxylase, which would not be obtainable from a 3 digit EC annotation (EC 6.3.4, Carbon—Nitrogen Ligases), which usually describes the chemical reaction.

In the future, a number of network diffusion limitations remain to be addressed. Here only enzymatic functions were considered, although ETA itself makes both enzymatic and non-enzymatic predictions using Gene Ontology (GO) terms [57]. The reason was that the network diffusion of labels taken from a GO directed acyclic graph (DAG) is more complex than from the simple EC hierarchy. Another concern is to further extend the coverage of yet unannotated proteins. As seen in Figure 5, ETA network diffusion performs better than a BLAST search when there are fewer homologs at high confidence z-scores. However, many non-homologous proteins share molecular function as a result of convergent evolution, [91] and variations can produce enzymes with similar function but differing sequence motifs [92]. Moreover, enzymatic function can be flexible and depend on context and expression level [7] such that enzymes are promiscuous and may perform several functions [93]. Presumably, to achieve even greater coverage, these problems will need to be addressed by raising the function detection sensitivity of the network. Further improvements in template construction [39], [94] or data integration [69] are possible directions towards these goals.

The nearest neighbor algorithm uses the same underlying network as the diffusion algorithm. It uses all ETA matches to proteins with known function and picks the function with the largest cumulative edge weights. The average edge weight of the winning function is used as a confidence value.

Coverage/accuracy curves were produced by sorting the predictions for the test sets in descending order by z-score, and then plotting a point showing cumulative accuracy (correct predictions/predictions made) and coverage (predictions made/size of test set) for each z-score threshold. The sensitivity/precision plot was produced in a similar manner: test set predictions were sorted in descending order by z-score and cumulative sensitivity (tp/(tp + fn) and precision (tp/(tp + fp)).

Comparison with BLAST was performed by BLASTing each chain in the test set against release 14 of Swiss-Prot, which is contemporary with the PDB data we used. The EC annotation of the resulting homolog with the smallest e value was taken as the predicted function. Comparison with PSI-BLAST was performed in an analogous manner against the same release of Swiss-Prot as the BLAST comparison. PSI-BLAST was allowed to run for 4 iterations.

The Structural Genomics testset was collected from the Protein Data Bank(PDB) by searching for proteins tagged with “structural genomics” and having an EC annotation. There were 1217 proteins with at least a 3-digit EC annotation, and 1036 with a 4 digit EC annotation. The list of candidate proteins for novel predictions were collected by searching the PDB for proteins labeled “structural genomics” and “unknown function.” Structures lacking annotation were counted by searching the protein databank for proteins tagged with “structural genomics” and matching a text search for “unknown function” or “hypothetical.”

The PDB 90 dataset was downloaded from the PDB website. Of approximately 18,600 proteins, 17,924 of them had enough homologs to perform an Evolutionary Trace. Annotations from the PDB were supplemented with annotations retrieved from the GOA database [95].

In the FLORA testset there were 911 domains which represent 821 unique PDB chains. Of those, Evolutionary Trace (ET) was able to find enough homologs for 806, which we used in our leave-one out experiments. The 15 proteins for which ET was not able to produce results are reflected in the lack of perfect sensitivity in Figure 2. ETA predictions of three-digit EC functions were made as described previously.

ClustalW [96] was used to calculate sequence identity between matches found between the testset of 1217 structural genomics proteins and their BLAST matches in SwissProt. Only matches with EC annotation were recorded. 10 sets of predictions were created, one for each {100, 90, 80, 70, 60, 50, 40, 30, 20, 10}% sequence identity cutoff. At each threshold we ignored matches that exceeded the allowed sequence identity. For BLAST predictions, the match with the lowest e-value below the given sequence identity threshold was used. Self matches were excluded.

For the S. aureus case study, the EFICAz² search was done through the PSiFR [97] tool provided at http://psifr.cssb.biology.gatech.edu/. The PSI-BLAST search was done through the web interface at http://www.ebi.ac.uk/Tools/psiblast/.

The structural comparisons were done with the DALI web interface [81]: http://ekhidna.biocenter.helsinki.fi/dali_server/.

Network distance calculations were performed via the networkX python library.

Table S1

Figure S1

Figure S2

Figure S3

Figure S4