Refine gene functional similarity network based on interaction networks

Background In recent years, biological interaction networks have become the basis of some essential study and achieved success in many applications. Some typical networks such as protein-protein interaction networks have already been investigated systematically. However, little work has been available for the construction of gene functional similarity networks so far. In this research, we will try to build a high reliable gene functional similarity network to promote its further application. Results Here, we propose a novel method to construct and refine the gene functional similarity network. It mainly contains three steps. First, we establish an integrated gene functional similarity networks based on different functional similarity calculation methods. Then, we construct a referenced gene-gene association network based on the protein-protein interaction networks. At last, we refine the spurious edges in the integrated gene functional similarity network with the help of the referenced gene-gene association network. Experiment results indicate that the refined gene functional similarity network (RGFSN) exhibits a scale-free, small world and modular architecture, with its degrees fit best to power law distribution. In addition, we conduct protein complex prediction experiment for human based on RGFSN and achieve an outstanding result, which implies it has high reliability and wide application significance. Conclusions Our efforts are insightful for constructing and refining gene functional similarity networks, which can be applied to build other high quality biological networks. Electronic supplementary material The online version of this article (10.1186/s12859-017-1969-1) contains supplementary material, which is available to authorized users.


Background
Most cellular components exert their functions through interactions with other cellular components [1]. The development of high-throughput measurement techniques such as tandem affinity purification, two-hybrid assays and mass spectrometry, has produced a large number of data, which is the foundation of biological networks [2]. Biological interaction networks, such protein-protein interaction network, gene regulatory networks, metabolic networks have been well studied and systematically investigated [3]. These networks play important roles in assembling molecular machines through mediating many essential cellular activities [4]. PPI networks occupy a central position in cellular systems biology and provide more opportunities in the exploration of protein functions in various organism [5,6].
In recent years, some researchers begin to pay their attention to the similarity networks, such as miRNA similarity networks [7][8][9][10], gene functional similarity networks [11,12]. Unlike the traditional interaction networks, similarity networks usually are constructed by measuring the similarity between the nodes in the networks. Since the similarity between each pair of nodes can be measured, these primary similarity networks usually are fully connected. For example, the construction of gene functional similarity networks is by measuring the sequence or ontology similarities between genes. The construction of miRNA functional similarity network is based on the functional similarity of two miR-NAs, which can be inferred indirectly by means of their target genes.
However, these fully connected similarity networks have one serious drawback. They do not meet the characteristics of biological network since they are fully connected [13]. Many previous studies have observed that biological networks are generally scale-free and their degree distributions follow the power law or the lognormal distribution [14][15][16]. From this point of view, we need to prune the unreasonable edges in the fully connected network. In the remainder of this section, we will first review some threshold selection methods, which have applied on gene functional similarity networks and phenotype similarity networks. Then we will put forward the proposed method.
Gene functional similarity networks have been widely used in some fundamental research, such as proteinprotein interaction prediction, disease gene identification and cellular localization prediction [11,[17][18][19]. Rui [11] constructed a gene functional similarity network to infer candidate disease genes on the genomic scale. The gene functional similarity network almost covers twice number of genes in the traditional PPI networks, which can enlarge the search range of candidate genes. However, the constructed gene functional network only keeps 100 nearest neighbors for each gene. As is pointed by Tian [20], this strategy is a very arbitrary for the selection of gene similarity values. Afterwards, Li [17] constructed a corresponding 5-NN network by means of keeping first five nearest neighbors of genes in the fully connected semantic similarity network. This method also has the common shortcomings with method Rui [11]. Besides, Elo [21] put forward a clustering coefficient-based threshold selection method to select a proper threshold for gene expression network. The similarity value below the selected threshold will be set to zero. However, small similarity in biological networks may be meaningful, while large similarity may also be noise. Perkins [22] applied the spectral graph theory on gene co-expression similarity networks for threshold selection. Perkins elaborated that applying a high-pass filter may remove some biologically significant relationships. These methods above always ignore the smaller similarity values, although they are meaningful sometimes.
At the same time, the threshold selection problem for the fully connected networks appears in other type of similarity networks [23][24][25][26]. For example, Van [23] made use of text mining method to classify over 5000 human phenotypes in the Online Mendelian Inheritance database and then constructed a fully connected phenotype similarity network. Li [24] employed the phenotype similarity network to infer phenotype-gene relationship. The authors only keep the first five nearest neighbors for each phenotype in the phenotype similarity network and obtain a 5-NN phenotype network. Later, Zhu et al. [25] come up with a new diffusion-based method to prioritize candidate disease genes. They believe that similarity values of phenotypes below the cutoff 0.3 are uninformative. Therefore, they did not considered similarity values below this selected threshold and set them to zero. Zou [27] and Vanunu [26] also keep the edge values higher than 0.3 in the phenotype similarity networks in their experiments. As for the phenotype similarity networks, the threshold selection has the same drawbacks with gene functional similarity network.
Based on the analysis for each method above, we can find that the threshold selection problem for the fully connected network is necessary, which has a significant effect on its applications. To the best of our knowledge, current threshold selection strategies for the fully connected networks are arbitrary or unreasonable. Therefore, it is still a challenge problem that how to construct a reliable gene functional similarity network.
In this article, we proposed a novel method to establish a high quality gene functional similarity network. The contribution of our study is listed as follow.
We construct an integrated gene functional similarity network based on six different functional similarity calculation methods. We built a referenced gene-gene association network based on the PPI networks. To the best of our knowledge, this is the first method that tries to refine gene functional similarity network based on a referenced genegene association network.

Methods
In this section, we will first introduce the experimental data briefly. Then we construct the integrated gene functional similarity network based on six functional similarity methods. After that, we will employ similarity indices between genes in PPI networks to construct nine gene similarity networks and get the referenced gene-gene association network. In the end, we obtain the refined gene functional similarity with the help of the referenced gene association network. Figure 1 depicts the flowchart of the proposed method.

Gene Ontology and Gene Ontology Annotation data
We downloaded the Gene Ontology (GO) data from the Gene Ontology database (dated July 2017) which contains 46,929 ontology terms totally subdivided into 4295 cellular components, 30,572 biological process and 12,062 molecular function terms. Gene Ontology Annotations (GOA) data for H. sapiens was downloaded from the Gene Ontology database (dated July 2017).

Protein-protein interaction data
Firstly, we obtain the protein-protein interaction data from human protein reference database (HPRD). HPRD is a high reliable PPI database, which is a resource for experimentally derived information about the human proteome. HPRD totally contains 39,240 interaction relationships relating 9617 proteins. Here, we select the maximum clique of HPRD, which contains 36,900 interaction relationships and 9219 proteins.
ConsensusPathDB are downloaded from the Website (http://consensuspathdb.org/). We selected three typical PPI networks based on ConsensusPathDB [28], which are Reactome, DIP and Biogrid. Specially, Biogrid contains 15,400 genes and 21,468 interactions, while Reactome contains 3332 genes and 19,604 interactions. As for DIP, it contains 3239 genes and 15,964 interactions. In this study, we will construct an integrated referenced gene-gene association network based on the four PPI networks above.

Construction of integrated gene functional similarity network based on GO and GOA
As we know, GO has three types of ontologies: cellular component (CC), molecular function (MF) and biological process (BP), respectively. Functional similarity between genes can be inferred from the semantic relationships of their annotated GOs [29]. Here we measure gene functional similarity using three types of ontology annotations that contain Inferred Electronic Annotations (IEA).
Since one method may have error prone in measuring functional similarity, the similarity here is calculated by six different kinds of methods. They are Resnik [30], Wang [31], GIC [32], SORA [33], WIS [34] and TopoIC-Sim [35] respectively. Method Resnik, Wang, and TopoICSim are pair-wise approaches, while method GIC, SORA and WIS are group-wise approaches. Besides, with the help of online tools [36,37], we can measure the gene functional similarity efficiently. In this article, 'functional similarity' refers to the similarity between genes, and 'semantic similarity' refers to the similarity between two GO terms.
Suppose there are genes A and B, the functional similarity between genes A and B can be measured from CC, MF and BP ontologies. Therefore, the functional similarity of gene A and Bis the integration of the three types of functional similarity, which can be measured by Eq. (1).
As for method Resnik, Wang, GIC, SORA, WIS and TopoICSim, their functional similarity results need to be integrated. The integrated functional similarity between genes A and B is calculated as follow: where MergedSim n (A, B)(n = 1, 2, 3,4,5,6) denotes the functional similarity method derived from method Resnik, Wang and GIC, SORA, WIS, TopoICSim, respectively.
Applying this operation to all gene pairs, thus we construct the integrated gene functional similarity network. It is noteworthy that the integrated gene functional similarity network is a fully connected network, which we need to purify the spurious edges in it. The number of genes in the integrated gene functional similarity network and PPI network is the same.

Construction of the referenced gene-gene association network
Here, we will construct a referenced gene-gene association network based on four PPI networks In order to maintain the unity of the number of genes, the genes in Reactome, DIP and Biogrid are the same with that in HPRD. We construct an integrated PPI network based on Reactome, DIP and Biogrid data in ConsensusPathDB and HPRD data. The construction process mainly has three steps.
Step one: construction of the weighted gene-gene association network We assess the reliability of protein-protein interactions in the integrated PPI network by edge clustering coefficient (ECC). Edge clustering coefficient is such a measure, which can both evaluate the reliability of interactions in PPI network and describe the association strength of two proteins [38]. For an edge E x, y connecting genes x and y, the ECC of edge E x, y is defined as.
where z x, y represents the number of triangles that actually include the edge in the network. d x and d y are the degrees of genes x and y, respectively. min(d x − 1, d y − 1) denotes the number of triangles that contains the edge E x, y at most. Obviously, the value of ECC(x, y) ranges from 0 to 1. Each pair of protein-coding genes in the integrated PPI network can be measured using Eq. (3), and we can obtain a weighted gene-gene association network.
Step two: construction of gene topological association networks For each pair of genes x and yin weighted gene-gene association network, a similarity score s xy is assigned to weigh their topological similarity. As we know, a higher similarity score corresponds to a higher probability of forming an association between two genes. Here, we define six similarity indices between two genes in the weighted gene-gene association network, which have been proposed by Yang [39]. They are the Weighted Common Neighbors (WCN), Weighted Resource Allocation (WRA) and Weighted Adamic-Adar (WAA) indices, as well as reliable-route weighted similarity indices [40,41]. The six similarity indices between genes x and y are formulated as follows: (1) Weighted Common Neighbors. where O xy denotes the common neighbor set of genes x andy, w xy represents the weight of the edge linking genes x and y, s z denotes the sum of weights for edges linking toz. Then, we will define another three similarity indices. Quasi-local similarity indices [42] not only consider the local similarity of two nodes, but also take local paths between them into account. Therefore, we define weighted reliable local path similarity indices as the similarity metric between unconnected genes x and y. The weighted reliable local path similarity indices are formulated as follows: where Γ(x) denotes the neighbor set of gene x, α is a parameter to adjust the contribution of length-3 paths. In this research, we set α as 0.5 to balance the length-3 path.
Applying those nine similarity indices to all gene pairs, we construct nine gene topological association networks, respectively. The edge values in the topological gene association networks denote the topological similarity between gene pairs.
Step three: construction of the referenced gene-gene association network By means of integrating the similarity scores in the nine gene topological association networks, we can obtain an integrated gene topological association network, whose edge weight is defined as.
where w i denotes the similarity score of gene pair in the ith gene topological association network. α i is the parameters to weight the nine gene topological association networks. α i was set as 1/9 to equally weigh the importance of the nine gene topological association networks.. In this article, we call this integrated gene topological association network as the referenced gene-gene association network. The edge values in the referenced genegene network denotes the topological similarities between gene pairs. The construction for the referenced gene-gene association network is completed.

Threshold selection for the integrated gene functional similarity network
Next, we will refine the integrated gene functional similarity network based on the referenced gene-gene association network. For any two genes A and B, their similarity values in integrated gene functional similarity network (IGFSN) and the referenced gene-gene association network (RGAN) are represented as sim(A, B) IGFSN and sim(A, B) RGAN , respectively. The similarity value between gene A and B in the refined gene functional similarity network (RGFSN) is denoted as sim(A, B) RGFSN , which can be calculated by Eq. (4).
Applying this operation to all gene pairs in the integrated gene functional similarity network, we can obtain the refined gene functional similarity network (RGFSN). From the Eq. (4), we can find that if the difference of similarity value between genes A and B in IGFSN and RGAN is large, the similarity value of A and B in RGFSN will be set to 0. In other words, the similarity value in IGFSN is noise according to RGAN. In this way, we can remove all the spurious edges in IGFSN.
What's more, taking the depth-first traversal experiment on RGFSN, we find that the refined gene functional similarity network have some isolated genes. The experiments results show that 8501 genes are formed one cluster, while the other genes (264) are isolated from this biggest connected component. As for this type of genes, we decide to add one of their neighbors in the integrated gene functional similarity network, to make RGFSN become one connected graph. At last, we can obtain a connected refined gene functional similarity network called RGFSN.
It is noteworthy that the small similarity value in integrated gene functional similarity network can be reserved based on our proposed method. Comparing with other threshold selection methods which filer out all edges with low similarity values, our method may be more reasonable.

Results
In this section, we will firstly compare the distributions of functional similarity values of different methods. Then we investigate the relationship between functional similarity values and protein proximity scores. After that, we focus on the global topological properties and the degree distribution of RGFSN. In the end, we conduct protein complex prediction experiment based on RGFSN, for verifying its reliability and application significance.

The distribution of functional similarity based on different methods
It is well accepted that gene functional similarity calculation methods used in this research have drawbacks [43]. For example, method Resnik has the 'shallow annotation' problem, while method Wang fixes the edge value of semantics contributions [31]. As for method GIC, it simply sums up the IC of terms when it measure the IC of a term set. Therefore, we propose a method to integrate the similarity results of the six methods to avoid the shortage of single method.
We investigate the distribution of six functional similarity methods and the integrated method. We randomly select ten hundred pairs of genes and then measure their functional similarity using method Resnik, TopoICSim Wang, GIC, SORA and WIS. The integrated functional similarity are computed by Eq.
(2). The distribution of functional similarity for the four methods are shown in Fig. 2.
From the results, we can clearly find that the highest functional similarity for method Resnik, GIC, WIS and SORA are not lager than 0.65, while the smallest similarity for method Wang is larger than 0.4. Obviously, this does not meet human perspective. By contrast, the integrated results are relatively reasonable. The highest and smallest functional similarity for integrated results are about 1.0 and 0.04, respectively. As a result, it is necessary for us to integrate the results of functional similarity methods.

Relationship of functional similarity and proximity scores
Next, we use the length of the shortest path between two genes in the integrated PPI network as their proximity measure. We choose 100 pairs of genes for each distance (one to five) and measure the functional similarity of gene pairs. To demonstrate the relationship between gene functional similarity scores and protein proximity scores, we draw the violin plot, which are shown in Fig. 3.
From the results, we can clearly find that gene pairs with closer distance (lower proximity scores) will have higher functional similarity scores. For example, the median functional similarity scores for distance one to five are 0.578, 0.519, 0.492, 0.475 and 0.458, respectively. The results indicates that the functional similarity scores are closely consistent with protein proximity scores. Therefore, we can construct a referenced gene-gene association network based on integrated PPI network to refine the gene functional similarity network. From this point of view, the proposed method is reasonable.

Global topological properties of RGFSN
The biological networks usually have their specific topological characteristics. We analysis the topology attributes of four networks based on Cytoscape 3.4 [44]. The corresponding results are presented in Table 1.
From the results, we can find that the topological properties of RGFSN meet the characteristics of biological networks, which are consistent with three other biological networks. For example, the diameter of a network refers to the longest distance between any two nodes [45]. The diameter of RGFSN is 8, while the diameters for HPRD, BioGRID and DIP networks are 14, 8 and 10, respectively. Besides, the cluster coefficient is a measure of the local interconnectedness of the network, whereas the path length is an indicator of its overall connectedness [46]. For biological networks, the cluster coefficient values are usually in the range 0.1 to 0.5 [47]. The cluster coefficients for HPRD, BioGRID, DIP, RGFSN are 0.102, 0.106, 0.098, and 0.118, respectively. Overall, RGFSN well meets the topological properties of biological networks.

Degree distribution of RGFSN
As is mentioned in previous section, many studies have observed that biological networks are generally scalefree. Their nodal degree distributions usually follow the power law or lognormal distribution [13,16] [48]. Here we employ four different models to fit the distributions of these four biological networks. These models are Gaussian distribution, power law distribution, log- Fig. 2 Distribution of functional similarity based on seven different methods. We can find that result for single gene functional similarity method is bias, while the similarity values for the integrated method are distributed from 0 to 1 evenly normal distribution and exponential distribution. All the fitting experiments are conducted on Origin 9. The results are shown in Table 2. Besides, the graphic view of the degree distributions for networks is shown in Fig. 4.
The detailed parameters (P) of four fitting models are listed in Table 2. The performances are evaluated by R-squares (R 2 ), which provides a measure of how well the data fits a certain model. The results show that RGFSN fits power law distribution best which is followed by exponential distribution. The R 2 scores for these two models are 0.9946 and 0.9816, respectively. As for BioGRID network, it fits the power law distribution best, while DIP and HPRD networks fit the exponential distribution best. From the results about the degree distributions, we can find that RGFSN has the typical characteristics of biological networks, e.g. scale-free, small world, rather than that of random network.

Protein complex detection experiment
Protein complexes are groups of associated polypeptide chains whose malfunctions play a vital role. Traditional methods predict protein complexes from protein-protein interaction networks, while some others are based on weighted association networks [43]. Here, we employ CPL [49] algorithm to predict protein complex based on RGFSN.
We verify the effectiveness and rationality of RGFSN by means of assessing the quality of predicted complex. To evaluate the clustering result, we used the jaccard score, which defined as follows: where K is a predicted cluster and R is a reference complex. Beside, we estimate the cumulative quality of the cluster result and set the MachScore as 0.25   [50]. Assume a set of reference complex R = {R 1 , R 2 , R 3 , ⋯, R n } and a set of predicted complex P = {P 1 , P 2 , P 3 , ⋯, P m }, the recall, the precision and F-measure at complex level are defined as follow.
A good prediction result should have higher accuracy, recall and F-measure values. The evaluation metrics about the quality of predicted complex have been  . 4 The graphic view of the degree distributions for each network discussed in detail [50,51]. In addition, the reference complexes was downloaded from CORUM database [52]. The number of reference complexes for human in this database is 1850 (see Additional file 1). We construct the 5NN network by keeping five nearest neighbors for each gene in IGFSN, which is proposed by Rui [11]. Here we call this network as the 5NN-IGFSN network. To increase contrast, we conduct the protein complex detection based 5NN-IGFSN with CPL algorithm. Besides, we also conduct protein complex prediction experiment based on HumanNet [53] and STRING [54] networks.
We evaluate the performance of CPL algorithm on STRING, HumanNet, 5NN-IGFSN and RGFSN according to the evaluation metrics. The results have been shown in Table 3. The precision, recall and F-measure of CPL algorithm based on RGFSN are 0.324, 0.347 and 0.314, respectively, while the results of precision, recall and F-measure for 5NN-IGFSN is 0.275, 0.223 and 0.246, respectively. From this point of view, the best performance in protein complex prediction indicates the reliability of RGFSN. The metric values for STRING and HumanNet are relatively low. The precision, recall and F-measure for STRING is 0.213, 0.268 and 0.236, respectively, while the results for HumanNet is 0.151, 0.142 and 0.146. Since many genes of HumanNet are not in CORUM database, its performance is worst.In the end, we take three examples to demonstrate the predicted results. Three referenced complexes are named as CNTF-CNTFR-gp130-LIFR, NCOR-HDAC3 complex and 20S proteasome, respectively. At the same time, we obtain three predicted complexes based on RGFSN using CPL algorithm. These three predicted complexes are shown in Fig. 5. The high overlap scores between prediction complexes and reference complexes demonstrate that RGFSN is a reliable biological network. The prediction results of CPL on RGFSN are presented (see Additional file 2).

Discussion and conclusions
In this study, we proposed a novel method to construct and refine the gene functional similarity network. Experimental results show that RGFSN is reasonable and effective. Thus, this method can be used to refine gene functional similarity networks effectively. However, two issues need to further study.

The construction of referenced gene association network
To refine the gene functional similarity network, we have to construct a reliable referenced gene-gene association network. This is the key point for the proposed method. In this study, we construct the PPI network that integrated four PPI data, which are DIP, Biogrid, Reactome and HPRD. The integrated PPI network is reliable and effective.
However, the integrated PPI network has itself shortcomings. It contains about 10,000 genes, which covers less than half of human genes. In addition, the integrated PPI network may be associated with false positives, although it has integrated many PPI networks. Therefore, we have to devote ourselves to seek other proper referenced network to achieve desired results in the next research.

The verification of the refined gene functional similarity network
How to verify the correctness and rationalization of RGFSN is a very challenging task. This is because there is no direct ways to evaluate the quality of the refined gene functional similarity network. In this research, we