Going the Distance for Protein Function Prediction: A New Distance Metric for Protein Interaction Networks

In protein-protein interaction (PPI) networks, functional similarity is often inferred based on the function of directly interacting proteins, or more generally, some notion of interaction network proximity among proteins in a local neighborhood. Prior methods typically measure proximity as the shortest-path distance in the network, but this has only a limited ability to capture fine-grained neighborhood distinctions, because most proteins are close to each other, and there are many ties in proximity. We introduce diffusion state distance (DSD), a new metric based on a graph diffusion property, designed to capture finer-grained distinctions in proximity for transfer of functional annotation in PPI networks. We present a tool that, when input a PPI network, will output the DSD distances between every pair of proteins. We show that replacing the shortest-path metric by DSD improves the performance of classical function prediction methods across the board.


Introduction
One of the best-studied classical problems on biological networks involves using proteins of known function, together with the structure of the network of known protein-protein interactions (PPI) to make predictions of functions of unlabeled protein nodes. This is an important problem because, even in the best-studied model organisms, such as S. cerevisiae, these networks contain many proteins whose function is still completely uncharacterized. There are many proposed methods for this problem, including versions of majority voting [1], neighborhood algorithms [2,3], clustering algorithms [4][5][6], algorithms based on maximum flow [7], or multi-way cut [8,9], and a number of others [10].
Modern approaches also seek to deal with data quality: generally, the known network is missing many valid interactions, and some interactions may be known with higher confidence than others [11,12]. Other modern approaches integrate PPI network data with high-throughput biological interaction data, such as sequence information, genetic interactions, structural information, and expression data [8,[13][14][15]. However, nearly all methods that predict function using PPI network structure depend, entirely or at least partially, on the simple shortest-path distance metric applied to the PPI graph.
We start with the observation that there is an intrinsic drawback to relying on the ordinary shortest-path distance metric in PPI networks. PPI networks are known to be ''small world'' networks in the sense that they are small-diameter, and most nodes are close to all other nodes (though the exact details of the degree distribution and the resulting extent to which they are ''scale free'' is a subject of lively debate, see [13,[16][17][18][19]). Thus any method that infers similarity based on proximity will find that a large fraction of the network is proximate to any typical node. In fact, this issue has already been termed the ''ties in proximity'' problem in the computational biology literature [4].
Furthermore, the fact that two nodes are adjacent (i.e., have shortest-path distance 1) in a PPI network can signify something very different than the adjacency of two other nodes. For example, as we discuss below, in PPI networks two nodes with many lowdegree neighbors in common should be thought of as ''more similar'' than nodes with few low-degree neighbors in common; and such nodes should also be thought of as ''more similar'' than two nodes whose common neighbors have high degree. Thus, characterizing node pairs based only on a shortest-path notion of distance fails to capture important knowledge encoded in the structure of the network.
What is needed instead is a finer-grained distance metric, capable of making more subtle distinctions of similarity than ordinary shortest-path distance would in a small-world network. To that end we introduce a new metric called Diffusion State Distance, or DSD. We show that DSD is much more effective than shortest-path distance for the problem of transferring functional labels across nodes in PPI networks.
We demonstrate the utility of the DSD metric by modifying a number of the most popular classical algorithms for function prediction to replace the ordinary notion of distance with DSD. For the problem of predicting functional labels in the S. cerevisiae network, this change improves all the algorithms we consider. We then show that similar improvements hold for the more sparsely annotated S. pombe network, implying that our advances should generalize to other biological networks.

Motivation for DSD
To start, we consider all known physical interactions in the S. cerevisiae PPI network -specifically, version 3.2.102 of BioGRID [20] on verified ORFs for S. cerevisiae, which contains 128,643 annotated physical interactions. After removing redundant edges and selecting the largest connected component, the resulting PPI network has 4990 nodes (where each ORF corresponds to a node) and 74,310 edges (where each edge corresponds to an annotated physical interaction). Figure 1(a) shows the histogram of shortest-path lengths from this network. The figure shows that almost all pairs of nodes (over 95%) are either 2 hops or 3 hops apart. Thus the ''typical'' distance between nodes under this metric is very small. Likewise, the overall network diameter is quite small as well.
The fact that most node pairs are quite close together means that the concept of ''neighborhood'' using this metric is not very useful. For example, if one looks at the 2-hop neighborhood of a typical node, it probably includes around half of all nodes in the graph. Hence, shortest-path distance cannot be expected to give a good sense of locality or modularity in this network, and we observe that if one seeks to use graph structure to infer similarity of function in a PPI network, using shortest-paths to measure distance has serious drawbacks.
One of the reasons that paths are typically short in biological networks like the PPI network is due to the presence of hubs -very high degree nodes which represent proteins that have many physical interaction partners. In fact, in the case of PPI networks, hubs often represent proteins with different functional roles than their neighbors. For example, chaperone proteins (proteins that help other proteins fold correctly) are often hubs, and are not   The correct functional annotation for GLR1, on the third level of the MIPS hierarchy, 32.01.01 (oxidative stress response) is found among none of its direct neighbors, but with the node that is closest in DSD, MXR2. MXR2 is closest in DSD because it has the most similar neighborhood to GLR1. doi:10.1371/journal.pone.0076339.g003 typically functionally related to their interaction partners. The same is true for proteins involved in translation. In our network, as an extreme example, the highest degree node is the protein NAB2, which has 2325 physical interactions in our PPI network. Its functional annotation terms: ''39-end processing'', ''RNA binding'' and ''RNA transport'' [21] suggest that this protein is involved in translation machinery and thus will bind with a highly diverse set of proteins with unrelated function. Hubs are also more likely to be proteins with multiple, distinct functions [22].
Hence, not all short paths provide equally strong evidence of similar function in PPI networks. Consider the network in Figure 1(b). Although nodes B and H are only one hop apart, this does not suggest they are functionally related, since H is a hub. Likewise, B and A are not necessarily functionally related, since they are connected through a hub.
To capture the notion that A and B are not necessarily related, we note that a random walk starting at A is not likely to reach B quickly. If we restrict attention to random walks of, say, 3 hops, then often one will not reach B from A at all.
This motivates the first element of our new metric definition. Given some fixed kw0, we define He fkg (A,B) to be the expected number of times that a random walk starting at A and proceeding for k steps, will visit B. Note that He fkg (A,B) bears some resemblance to previous diffusion approaches suggested for PPI networks, see [23,24], though we will next be building metric structure on top of He fkg (A,B) in a completely novel way. For now, note that He fkg (A,B) captures our intuition regarding similarity, because node pairs connected by many short paths of low-degree nodes will tend to have high He fkg () values. If node pairs with a large He fkg () value are then somehow considered 'similar' in our metric, then clearly He fkg () does a better job of capturing 'similarity' than does shortest-path distance. This is because A and B are relatively far under this metric; the influence of the hub H decreases the likelihood of a random walk from A landing at B.
But a metric whose notion of similarity is based only on using He fkg () between the two given nodes directly does not solve all our problems. In particular, note that He fkg (A,H) will indicate strong similarity, even though (as we have argued) A is not likely to be strongly functionally related to H. Furthermore, He fkg () is still far from a metric; note it is not even symmetric. Also, He fkg () is only a pairwise, and not a global measure of node similarity. This observation leads us to use He fkg () in a more subtle way, and in a different manner than previous diffusion approaches, resulting in the definition of DSD, which we describe next.

Definition of the New Distance Metric
Consider the undirected graph G(V ,E) on the vertex set V~fv 1 ,v 2 ,v 3 ,:::,v n g and jV j~n. Recall that He fkg (A,B) is defined as the expected number of times that a random walk starting at node A and proceeding for k steps, will visit node B. In what follows, assume k is fixed, and when there is no ambiguity, in the value of k, we will denote He fkg (A,B) by He (A,B). We further define a n-dimensional vector He where jjHe(u){He(v)jj 1 denotes the L 1 norm of the He vectors of u and v.
We will show for any fixed k, that DSD is a metric, namely that it is symmetric, positive definite, and non-zero whenever u=v, and it obeys the triangle inequality. Thus, one can use DSD to reason about distances in a network in a sound manner. Further, we show that DSD converges as the k in He fkg (A,B) goes to infinity, allowing us to define DSD independent from the value k. Figure 2 shows the distribution of DSD values in the PPI network of S. cerevisiae as downloaded in BioGRID [20]. The figure shows that DSD distances take on a smooth, continuous range of values. This is in sharp contrast to shortest-path distances, as shown in Figure 1(a), and shows that DSD is capable of making much finer-grained distinctions of similarity over pairs of nodes in the network.  GLR1's immediate neighbors contain its correct functional label at the third level of the MIPS hierarchy (32.01.01: oxidative stress response). However, the node closest in DSD is MXR2, which has exactly the same functional label. DSD recognizes that nodes having large common low-degree neighborhoods are highly similar and correctly identifies functionally similar node pairs, and does so in situations where shortest-path distance fails.

Functional Categories
We continue to work with the dataset described above, the largest connected component from an experimentally derived physical interaction PPI network from S. cerevisiae with 4990 nodes, and 74,310 edges. As in the papers introducing the classical function prediction methods we consider, we primarily use the MIPS (Munich Information Center For Protein Sequences) functional catalogue (FunCat) for our functional category labels [21]. We use the latest version of FunCat (version 2.1) and the first, second and third level functional categories, retaining those for which at least three proteins in our dataset are annotated with those labels. We present results for the first level (17 functional categories) second level (74 of the 80 functional categories that annotate at least one protein, annotate at least 3 proteins) and third level (154 of 181 functional categories that annotate at least one protein, annotate at least 3 proteins) MIPS annotations. MIPS is a shallow, leveled, hierarchical classification scheme for protein function, but we also present results for the popular Gene Ontology (GO) [25], where the variable depth hierarchy of the annotation labels makes the evaluation of labeling methods more complicated.
We assumed all published labels were correct and did not attempt for this study to weigh them by a level of confidence in each type of biological experiment. The following classical methods were tested in their original form (using shortest-path distance), against a DSD variant in cross validation.

Cross Validation Task
We considered 2-fold cross validation tasks. In each of the 2-fold cross validation tasks, we first randomly split the annotated proteins into 2 sets, and consider only the annotations on one partition in the training set, when trying to predict the annotations on proteins in the test set, and average the performance over the 2 folds of the cross validation. We conduct 10 runs of 2-fold crossvalidation and report the mean and standard deviation of the following performance measures over these 10 runs.
Accuracy. This is the performance measurement suggested in [1]. Each protein is assigned its highest-scoring functional label. The label is considered correct if it appears among the labels that were assigned to the protein. We calculate the percentage of proteins that are assigned a correct label.
F1 score. This is the performance measurement suggested in [26]. For each protein, the possible functional labels the algorithm could assign are stored in a ranked list according to score. Each label is considered correct if it appears among the labels that were assigned to the protein, and incorrect otherwise. We calculate precision and recall by looking at the top a (in our case, we present results for a~3) predicted annotations. Then the F1 score for each function can be calculated as: We average F1 scores over the individual functions and obtain the overall F1 score for each algorithm.

Protein Function Prediction Methods
Neighborhood Majority Voting Algorithm. This is the simplest of all function prediction methods. Directly applying the concept of 'guilt by association', Schwikowski et al. [1] considered for each protein u[V its neighboring proteins. Each neighbor  votes for their own annotations, and the majority is used as the predicted functional label. To incorporate DSD, the neighborhood of u is defined simply as the t nearest neighbors of u under the DSD metric. Furthermore, two schemes are considered: an unweighted scheme where all new neighbors vote equally, and a DSD weighted scheme where all new neighbors get a vote proportional to the reciprocal of their DSD distance.
x 2 Neighborhood Algorithm. In the original x 2 neighborhood algorithm [3], each annotation a present in protein u's neighborhood will be assigned a x 2 value based on the following formula: where n a is the observed number of annotations a in the neighborhood, and e a is the expected number of annotations a based on the whole protein-protein interaction network. Then protein u is assigned the functional label a with the maximum x 2 value. Again, it is straightforward to adapt this to use DSD: the neighborhood of u is simply defined as the t closest nodes to u under the DSD metric.
Multi-way cut Algorithm. We consider the minimal multiway k-cut algorithm of Vazquez et al. [9] as implemented by [7]. The motivation is to minimize the number of times that annotations associated with neighboring proteins differ. In particular, the dual ILP (integer linear programming) is formulated, so we instead seek to maximize where the edge variables X u,v,a are defined for each function a whenever there exists an edge between proteins u and v. It is set to 1, if protein u and v both are assigned function a, and 0 otherwise. The node variables X u,a are set to 1 when u is labeled with function a and 0 otherwise. The first constraint insures that each protein is only given one annotation. The second constraint makes sure only annotations that appear among the vertices can be assigned to the edges. While this problem is NP-hard, the ILP is tractable in practice; in our case we use the IBM CPLEX solver (version 12.4, dated 12/2011, http:// www.ilog.com/products/cplex/). For the DSD version, we simply add additional edges between vertices whose DSD is below a threshold. We set a global threshold D based on the average DSD of all pairs, specifically we set D~m{c Ã s, where m is the average, and s is the standard deviation of the global set of DSD values among all pairs of nodes in the graph. We experiment with c in the range f1:5,2:0,2:5,3g.
Functional Flow Algorithm. Nabieva at al. [7] use a network flow algorithm on the graph of protein interactions to label proteins. The idea is to consider each protein having a known function annotation as a 'reservoir' of that function, and to simulate flow of functional association through the network to make predictions. We adapt the approach to use DSD by creating an edge between each node pair, with a weight inversely proportional to DSD. For computational efficiency we do not create edges when the reciprocal of DSD is below a small value. As in the original functional flow, we calculate flow through this new network at each time step. We denote the size of the reservoir of function a at node u and time step i, to be R a i (u). For a given function (annotation) a, we initialize the reservoir size at node u to be infinite if protein u has been annotated with function a; otherwise we set it to be 0. More formally: , if u is annotated with a 0 otherwise We then update the reservoir over a sequence of timesteps (we use 6 timesteps, as in the original version): where g a t (v,u) is the amount of flow a that moves from u to v at time t. We incorporate DSD into the edge weight as follows:  The final functional score for node u and function a over 6 timesteps is computed as the total amount of incoming flow.

Formal Properties of DSD
We now present the formal proofs that DSD is a metric for any fixed k, that it converges as k goes to infinity, and obtain the explicit form that allows the calculation of DSD values in the limit. Lemma 1. DSD is a metric on V, where V is the vertex set of a simple connected graph G(V ,E).
Proof. Clearly the DSD of a node to itself is 0, and DSD is non-negative and symmetric. It remains only to show that DSD satisfies the triangle inequality, and that DSD(u,v) is strictly positive whenever u=v.
For all u,v,w[V , we have by the L 1 norm property: {He(w)jj 1 and therefore: Thus, the triangle inequality follows easily: Next we prove the identity of indiscernibles, namely, DSD(u,v) is non-zero for all u=v. We first need some notation. We define the one step transition probability matrix P as the n-dimensional square matrix, whose (i,j)th entry is given by: where d i is the degree of node v i Vi~1,2:::,n. Note that P represents the probability to reach each neighbor in the random walk where all neighbors are reached with equal probability. Then the definition of k step transition probability matrix P fkg~Pk follows for all positive k. We also define the 0 step transition matrix P f0g to be the identity matrix I, because every random walk starts from the source vertex, and thus a zero-step random walk reaches its source vertex with probability 1 and all other vertices with probability 0. We denote by p flg ij or p flg vi,vj the (i,j)th entry in the lth transition probability matrix, where l §0 and i,j[f1,2,:::,ng.
For each ordered pair of vertices (v i ,v j ), recall that He(v i ,v j ) is defined as the expected number of times that a random walk with length k starting from v i will visit v j . In order to calculate He(v i ,v j ), we define an indicator variable I Now we are ready to show the identity of indiscernibles. Recall we are assuming the graph is connected. It is trivial when n~jV j~2, so consider the case where n §3. We prove this next by contradition.
Assume Since the zero-step transition matrix P f0g is an identity matrix, we know that p f0g vi,vj~1 for i~j and 0, otherwise.
Thus we have where the third line represents a set of n{2 equations: one for each v i , distinct from a and b.
By multiplying a factor p vi,a~p f1g vi,a to both sides of all equations in the third line and summing over i, we have: By completing the sum over i, we have: By applying to the equation above the first and the second equation in (1), where: we have: [½0,1, we have p fkz1g a,a~1 . We next argue that it must be the case that d(v i )~1, for all v i with (v i ,a)[E, namely, all of a's neighbors must have degree 1.
Starting from , and Vv i : (v i ,a)[E and p fkg a,vi w0, p vi,a~1 . Since we have X (vi,a)[E p fkg a,vi~1 , there must exist at least one neighbor x, s.t. p fkg a,x w0, and p x,a~1 as well. Because p x,a~1 , we know that a is the only neighbor to x and d(x)~1. Also due to the fact that p fkg a,x w0, there exists a path of length k from a to x, which should consists of two parts, one path (denoted as path Ã (a?a)) of length k{1 from a to a and an edge from a to x as the last step in the path.
Consider any neighbor y : (y,a)[E. By using the path Ã (a?a) and the edge (a,y)[E, we can construct a path of length k from a to y and therefore p fkg a,y w0. Thus, p y,a~1 , and therefore d(y)~1. Thus we have shown that it must be the case that d( As a result, all of a's neighbors are only connected to a. If (a,b)[E, then a and b can't have any more neighbors because they must be the only neighbor of each other; thus for all c[V {fa,bg, c is not connected to a or b, which contradicts the fact that the graph G is connected. If (a,b) 6 [E, a and b are not connected because all of their neighbors are only adjacent to themselves respectively. Therefore, the existence of such pair (a,b) where a=b and DSD fkg (a,b)~0 contradicts the fact that the graph is connected, and we can conclude that DSD(a,b)~0 if and only if a~b.
In the above we were reasoning about DSD for a fixed value of k. Denote by DSD fkg the value of DSD for a particular fixed k. Next we discuss how DSD values depend on k. We show that when the one-step transition matrix P is diagonalizable then DSD fkg (u,v) converges to a stationary value DSD(u,v).
Lemma 2. Let G be a connected graph whose random walk one-step transition probability matrix P is diagonalizable and ergodic as a Markov chain, then for any u,v[V ,DSD fkg (u,v) converges as k, the length of the random walk, approaches infinity.
Proof. Since the one step transition probability matrix P is diagonalizable, we have by diagonalization: where the orthogonal matrix V~U {1 ; each column of V is the normalized right eigenvector of P, each row of U is the normalized left eigenvector of P, and L~diag(l 1 ,l 2 ,:::,l n ) is the diagonal matrix of all eigenvalues where 1~l 1 wjl 2 j §::: § jl n j by ergodicity. Let v 1 ,v 2 ,:::,v n denote the normalized orthogonal column vectors in V and u T 1 ,u T 2 ,:::,u T n the normalized orthogonal row vectors in U. We denote v ij as the jth entry in vector v i and u ij as the jth entry in vector u i .
Thus it follows that P ftg~V L t U, namely Vi,j[V ,t §1: Next, for all i,j,b[V , we define an infinite sequence A t (i,j,b) for t~1,2,:::,?: where I is the identity matrix. Therefore, by definition, we can rewrite DSD fkg as the partial sum: Before we prove the claim, we show that it produces what we need. If the claim is true, then the limit j lim k?z? X fkg t~0 A t (u,v,b)j exists, and we can denote it as A Ã (u,v,b). Then it follows that: and we have proved convergence.
But Claim 1 is true from the Cauchy root test, namely,  (4) holds because 1~l 1 wjl 2 j §::: § jl n j and jjv c jj 2~j ju c jj 2~1 . We now use this lemma to produce an explicit way to calculate in the limit. This frees DSD entirely from dependence on the parameter k, and this is the version of DSD we use in our experiments. (We remark, that for our yeast PPI network, we found empirically that DSD fkg values were very close to the limit when k §5.) Lemma 3. Let G be a connected graph whose random walk one-step transition probability matrix P is diagonalizable and ergodic as a Markov chain, then for any u,v[V , we have , where I is the identity matrix, W is the constant matrix in which each row is a copy of p T , p T is the unique steady state distribution, and for any i[V , b T i is the i th basis vector, i.e., the row vector of all zeros except for a 1 in the i th position.
Proof. To start, we denote the matrix P{W by D. By Lemma 2 where we have shown that lim k?? DSD fkg (u,v) exists and is finite, we can denote the limit as simply DSD(u,v) in the following context. It follows that for any u,v[V : E lim k??
where (6) holds by definition, (7) holds because b T u W~b T v W~p T , (8) holds because (P n {W )~(P{W ) n for all n §1 by Claim 2, and (9) holds by Claim 3, where we finish the proof of Lemma 3 by proving Claims 2 and 3 next.
Claim 2. (P n {W )~(P{W ) n , for any non-zero positive integer n.
Proof. The proof is based on the following observations: WP~W and PW~W , so WP k~W and P k W~W . Further, W m~W , so W k P m~W and P m W k~W for any integers kw0 and mw0. Now, if you construct the binomial expansion of (P{W ) n , each cross term of the type W k P m or P m W k reduces to W , and all of the W s cancel out except for one, leaving (P n {W ).
exists. We show that (I{D) has full rank by showing that 0 is the only vector in the left nullspace of (I{D). That is, if x T (I{D)~0, then x~0.
x T (I{D)~0 x T {x T D~0 x T D~x T x T D k~xT for anytw0 Now, we already know that lim k?? D k~l im k?? P k {W exists and is 0 by ergodicity. So we have lim t??
x T D k~xT x T 0~x T which is only satisfied by by x~0.
1. Since (I{D) {1 exists and lim k?? D k~0 , we can calculate as follows: where C is defined as Iz lim k?? X k t~1 D t , which exists because (I{D) {1 exists and lim k?? D k~0 .
That finishes the proof for Lemma 3.

MIPS Results
Both the DSD majority voting method and the x 2 neighborhood method sort vertices in order of their smallest DSD to a given vertex, and set a threshold t, where the first t vertices participate in the functional vote. Table 1 gives the results in 2-fold cross validation for both these methods when we set t~10, whereas Figure 4 gives more details about the dependence on t; basically once enough neighbors were included (i.e. t §8) results of the DSD version of majority voting converged. Similar details about how the x 2 neighborhood method depends on t appears in ( Figure S1). For the multi-way cut and functional flow results, we considered vertices to be in the DSD neighborhood of a node v if their DSD from v was less than c standard deviations below the mean DSD value across the entire dataset. Table 1 presents the 2-fold cross validation DSD multi-way cut results and DSD functional flow results with c~1:5. In addition, Figure 5 gives more detail about the dependence of DSD functional flow on the parameter c, where we tested c[f0:5,1:0,1:5,2:0,2:5,3g. Similar details about how the multiway cut method depends on c appear in ( Figure S2).
We notice from Table 1 that in every case, the DSD versions of all four methods perform better than the versions based on ordinary distance. It is particularly interesting that using DSD improves the functional flow algorithm, since functional flow already seeks to capture a notion of diffusion in the graph. In fact, the DSD versions of Majority Voting perform better than all four of the methods based on ordinary distance. The best performing algorithms among all the ordinary distance algorithms on the yeast MIPS annotations at all three levels of the MIPS hierarchy were the majority voting algorithms. The best performing algorithms overall are the DSD version of the majority voting algorithms, and they clearly achieve the best performance compared to all four ordinary distance and all three other DSD-based methods. In particular, they achieve 63.7% mean accuracy and 47.2% mean F1 score (unweighted) and 63.2% mean accuracy and 48.1% mean F1 score (weighted) on the first level of the MIPS hierarchy. This is compared to the original majority voting algorithm, which only obtains 50.0% mean accuracy and 41.6% mean F1 score. Hence DSD provides more than 13% improvement in accuracy and more than 5% improvement in F1 score for MIPS top level functional categories, and improves on ordinary DSD on the second and third levels of the MIPS hierarchy across the board.

GO Results
GO (Gene Ontology) [25] is a deeply hierarchical classification scheme, which makes defining and evaluating function prediction methods much more complicated. For this reason, most function prediction methods for yeast PPI networks use the ''flatter'' set of MIPS categories, but GO is the more widely used ontology. Thus we wished to measure performance for GO functional labels, despite the difficulties. For example, in GO, sometimes, nodes are labeled with child functions but not labeled with their parent functions, even though it is assumed that child functions are more specific and inherit the functions of their parent nodes. How much credit should be given to a node when we do not label it correctly at its most specific level, but succeed in labeling it with a less specific ancestor term? The exact match and functional path methods of Deng et al [27,28] are designed to directly answer this question and evaluate methods that perform GO annotation. We tested the performance of ordinary distance versus DSD versions all four protein function prediction methods considered above also on the GO, and evaluated the results using the exact match method and functional path method of Deng et al. [27]. We find, again, that using the DSD-based algorithm improves performance.
We consider labels in the biological process category of the GO hierarchy (we used OBO v1.2 [29], data version 2013-07-18) along with annotations downloaded from the SGD database (data version 2013-07-13). We exclude GO terms that are annotated with evidence codes ''IEA'' ''RCA'' or ''IPI''. For each protein that was labeled with a term in the GO hierarchy, we automatically also label it with all more general parent terms in the GO hierarchy as well. We then, define the informative nodes in the GO hierarchy to be functional annotation terms that 1) are at least three levels below the root and 2) are terms that annotate more than 50 proteins in our dataset, where the second condition on informative nodes, and the number 50 is suggested by the method of Deng et al. [27]. It is these GO functional terms, somehow capturing the middle levels in functional specificity, that we use instead of a level of the MIPS hierarchy for the functional labels. The result is 136 informative biological process GO terms among 4322 out of 4990 ORFs that will comprise our labeling set (the total number of annotations is 58,519); thus we have a label set that is close in size to the one from the third level of the MIPS hierarchy. The ''exact match'' method then simply counts the number of correct labels, just as we did for a fixed level of the MIPS hierarchy. We see similar improvements for GO annotation as we did for MIPS annotation; detailed results for accuracy and F1 score for all four methods appear in ( Figure S3).
However, note that this exact match evaluation gives no credit when we label a protein with an incorrect child that still has many ancestor terms in common with the correct label. Thus, for a more fair evaluation of the predictions that takes into account hierarchical relations among the GO labels, we also use the functional path method of Deng et al. [27]. This presents a plausible way to give partial credit when a node is labeled partially correctly, in the sense that the lowest depth child label assigned to the node and the correct label of the node are different, but if we consider the path in the ontology from the root to these two labels, there are ancestor labels these two functional paths have in common. We use exactly the functional path method of [27] to calculate precision and recall values for each protein, and then overall precision and recall values are averaged over all the proteins. We calculated both alternative ways to count functional paths overlap presented in Deng et al. [27], one that simply counts the number of nodes that appear jointly in sets of known and predicted functional annotations (which we will call the overlap counting method) and the other which takes into account at what depth the functional paths from both sets diverge (which we will call the overlap depth method).
Setting t, the neighborhood threshold for DSD majority voting anywhere in the range starting from 2, the DSD version of majority voting improves precision and recall simultaneously, as compared to original majority voting, regardless of the method used for counting overlaps. For example at t~8, we get precision = 67.3% 60.3% and recall = 40.3% 60.2% compared to precision = 59.6% 60.4% and recall = 34.3% 60.4% for the overlap counting method, and we get precision = 61.1% 60.3% and recall = 27.3% 60.1% compared to precision = 52.9% 60.5% and recall = 20.5% 60.4% for the overlap depth method, for 2-fold cross validation on the yeast PPI network. In fact, we find that all DSD versions of the four algorithms perform better than their ordinary distance versions, regardless of the methods used for keeping track of overlaps. Figure 6 shows the improvements in F1 scores (setting t~10 and c~1:5 as before) over all three methods of counting the overlaps. Thus DSD is improving functional annotation also for GO annotation.

Discussion
Our function prediction results demonstrate the utility of defining and using a fine-grained distance measure that is specifically tailored to the subject application domain. For the PPI network, the observation that shortest-paths that go through hubs are less informative led to the particular design of the He () measure incorporated into the DSD definition. The He () measure has connections to other diffusion-based measures previously proposed [23]; however, using He () globally across all nodes, and comparing such vectors via L1 norm is entirely new and enables DSD to capture truly global properties of network topologies.
We have shown that substituting DSD for ordinary shortest path distance results in dramatic improvements when using network information to predict protein function for the S. cerevisiae PPI network. We expected that similar improvements would hold for the PPI-networks of other organisms as well, though the sparseness of current experimental known functional annotation and the number of PPI interactions currently known for other organisms means that we would not yet expect absolute levels of accuracy that are comparable to those achieved for S. cervesiae. We tested this intuition by also considering what is known of the PPI network of a different, less-well annotated yeast species, S. pombe. Based on the interactions in BIOGRID [20] version 3.2.102, the largest connected component has 1925 nodes and 4874 edges. Thus compared to the S. cerevisiae network, a much smaller and sparser subset of the PPI network is known. Figure 7 shows the shortest path and DSD distance distribution of this network. As expected, it is low-diameter but not yet as low-diameter as the S. cerevisiae network, and the DSD distribution is more spread out, but not yet as smooth as for the S. cerevisiae network. We would predict that the distribution would start looking more like S. cerevisiae as more of the network becomes known. We found no reliable MIPS annotation for S. pombe, so we instead used GO annotation, which we downloaded from the Pombase database, data version 2013-07-01 (http://www.pombase.org). We extract ''biological process'' GO terms only and remove GO terms annotated with evidence codes ''IEA'', ''IPI'' and ''RCA'', leaving 85 informative (see GO results section above, for definition) GO terms annotating 1722 out of 1925 proteins, with a total number of 5818 annotations. We did the ''exact match'' version of GO evaluation for 10 runs of cross-validation, and looked at the difference in performance for majority vote and functional flow methods using ordinary shortest path distance and DSD. We see similar improvements using DSD as with S. cervisiae; results appear in ( Figure S4).
We present a simple to use, freely-available tool that given any PPI network will produce the DSD values between each pair of proteins, either as a webserver or for download from http://dsd.cs. tufts.edu/. This tool both allows the calculation of DSD fkg for some chosen k, and, as in the results presented in this paper, for DSD in the limit. In fact, we also tried DSD f5g in place of DSD for all the methods in this paper and results were already quite similar to DSD in the limit. We suggest based on the dramatic improvements for the S. cerevisiae PPI network that DSD values be used in place of ordinary distance when using network information to predict protein function, either using the PPI network alone, or as part of a modern integrative function prediction method that includes data from a variety of sources beyond the PPI network, such as sequence information, genetic interaction, structural information or expression data [8,[13][14][15].

Supporting Information
Figure S1 Improvement of mean Accuracy and F1 Score for DSD at different neighborhood thresholds for the x 2 neighborhood algorithm in 10-runs of 2-fold cross validation (with standard deviation error bars). (TIF) Figure S2 Improvement of mean Accuracy and F1 Score for DSD at different neighborhood thresholds for the multi-way cut algorithm in 10-runs of 2-fold cross validation (with standard deviation error bars). (TIF) Figure S3 Improvement of mean Accuracy and F1 Score for DSD at different neighborhood thresholds for all four methods in 10-runs of 2-fold cross validation (with standard deviation error bars) using the GO catagories (using the method that gives wcredit only for exact matches to each GO term). (TIF) Figure S4 Performance evalutaion for GO term prediction on the S. pombe network. (a) Comparison of mean F1 score of DSD and non-DSD majority vote (setting t = 10) and functional flow (setting c = 1.5) algorithms using all three methods of counting GO term matches; (b1,b2) Comparison of mean accuracy and F1 Score under the exact match method for different settings of the parameters t and c, over 10 runs of 2-fold cross validation (with standard deviation error bars).