FSM: A Fast Similarity Measurement for Gene Regulatory Networks via Genes' Influence Power

1 Introduction

Gene regulatory network (GRN) analysis has attracted increasing attention in bioinformatics and data mining, as some complex diseases like cancer or diabetes are usually caused by dysfunction of relevant networks or network communities rather than mutations of individual molecules [de Souza Jacomini et al., 2016]. GRNs illustrate the interactions between molecular regulators in organisms, such as DNA, RNA, proteins, and other chemicals. Researchers often use data mining techniques on GRNs to uncover the internal relationships as well as essential pathways of diseases, particularly cancer[Ruan et al., 2015; Denitto et al., 2015].

∗Corresponding Author

A central problem in GRNs is how to measure the similarities between two networks, each of which is under a specific condition. For example, given several GRN obtained from clinical data, some of which are labeled as specific cancers, can we have a "ruler" to measure the distance between any two of them and determine if the rest of the samples are also suffered from the same cancer or not (such as the example illustrated in Figure 1)? Or given a sequence of GRNs obtained from patients to be diagnosed, can we quickly scan them, capture the similarities and differences between them, and detect the anomalies? Moreover, by searching and measuring such similarities, the function of unknown molecules or biological pathways associated with specific diseases can be determined. Similarity problem in GRNs is more complex than it is in other domains such as social networks. In social networks, we focus more on the structural similarities and consider nearly no additional assumptions or properties on nodes and edges. However, when we consider similarity measurement problems in GRNs, we actually consider the functional similarity of two networks. Properties of nodes and edges are highly respected, including the properties which are handdesigned or obtained from biological experiments.

Various computational methods have been proposed to compare graphs by learning a distance function [Liu et al., 2019]. One widely used method is graph edit distance (GED). GED calculates the minimal cost of transforming a graph into another. Generally, optimization of the cost function is NPhard [Zeng et al., 2009]. another state-of-art method is graph kernel. Graph kernels compare the substructure of two networks by the following steps: first it divides networks into sets of characteristic substructure patterns, and the cut-set of the pattern sets is determined [Horv´ath et al., 2004]. Then, the similarity of two networks is calculated by kernel function k : ϕ(X) × ϕ(X′ ) → R. ϕ is a decomposition algorithm for graph X → H into a Hilbert space, such that k(X,X′ ) = I(ϕ(X), ϕ(X′ )), where I is a isomorphism algorithm for substructures. Graph kernels usually have good performance in some application of real-world networks, because they consider local or global properties of networks. However, finding such a similarity using graph kernels is usually a NP-complete problem [G¨artner et al., 2003].

Accordingly, computational complexity is a major challenge in graph similarity measurement. Existing methods are not usually suitable for practical applications in biology and clinical research because real-world GRNs are highly complex systems with thousands of nodes. To measure the similarity between GRNs faster and credibly, a method called fast similarity measurement (FSM) is proposed. This method is inspired by the concept of generalized hamming distance (GHD) which is capable of identifying subtle variations in the topology of paired graphs [Mall et al., 2017]. Besides, the concept of influence power (IP) is introduced to identify influential nodes in a GRN. FSM can fast measure the similarity of paired GRNs by detecting the differences in nodes' IP between two GRNs. The contributions of this study are as follows:

• FSM captures the subtle variations between influential nodes and provide a more credible GRN similarity measurement.
• The concept of nodes' influence power (IP) is introduced. It is inspired by real-world phenomena and can effectively recognize influential nodes among the GRN.

The rest of this paper is organized as follows. Section 2 introduces some preliminary concepts and describes IP and GHD. In Section 3, IP and the FSMmethod are described in detail. In Section 4, the method is evaluated on a large number of datasets of varying size and sources. Finally, Section 5 concludes the paper.

2 Preliminaries

2.1 Definitions and Problem Statement

Let G = (V, E) be a directed unweighted graph, where V is the node set and E is the edge set. In a GRN, each gene is represented as a node in V . And a directed edge between genes, e.g., (u, v) ∈ E implies that gene u regulates gene v in the regulatory network. For the node v ∈ V , inN k (v) is defined as the k-th order incoming neighborhood set of v (k = 1, 2, . . . ). When k = 1, for each node in the set there exists a directed edge pointing at v, i.e., inN 1 (v) = {∀u|(u, v) ∈ E}. Furthermore, inN 2 (v) = {∀u|∃w, (u,w) ∈ E, (w, v) ∈ E}. Similarly, the k-th order outgoing neighborhood set of v is denoted by outN k (v). That is, for all u ∈ outN k (v), the distance from v to u is k in the directed graph. Finally, |·| denotes the number of elements in a set. For example, the number of nodes in a GRN is |V |, and the in-degree and the out-degree of node v are ∣∣inN 1 (v)∣∣ and ∣∣outN 1 (v)∣∣, respectively.

Based on the definitions above, we formalize the problem as follows.

Problem 1. supposing GA = (V,EA) and GB = (V,EB) are two different GRNs sharing the same node set V , and the permutation of nodes keeps unchanged. The distance between the two networks is D(GA,GB), where D is a distance metric. Let p = P(y = 1|D(GA,GB)) be the conditional probability that GA and GB belong to the same class given the distance D(GA,GB), and (1 - p) = P(y = 0|D(GA,GB)) be the probability of the opposite situation. We aim to find a proper distance metric D, satisfying:

In our experiment, the dataset contains networks labelled as different classes (e.g. class1, class 2, ...). First, We randomly chose two networks with same scale from two different classes respectively, denoted as GA and GB. Then, a series of permutated networks (namely GA1,GA2, · · · and GB1,GB2, . . . )are generated by randomly swap some links in GA and GB. D(GAi, GA) and D(GBi, GA) (i = 1, 2, . . . ) are calculated. For the former case, we label y = 1, and for the latter case, we label y = 0. Then, a 10-fold crossvalidation SVM was performed to check whether D can correctly measures the similarity between the networks and satisfy the minimum entropy. The above process is repeated several times to ensure that as many samples as possible in the dataset are covered.

2.2 Influence Power and Generalized Hamming Distance

Since many real-world GRNs are thought to be scale-free [Caldarelli et al., 2002], a large proportion of genes are regulated by only a few hub genes. If each gene is thought of as a person in a community, and the regulatory relationship between genes as exchange of economic goods, the differences in income and expenses will result in social class inequality. A gene's "social class" is associated with the level at which it regulates other genes and the level at which it is regulated by other genes. Here, a metric called influence power (IP) is proposed to quantify the degree of regulation based on the topological structure of GRNs. IP measures a node's influence in terms of its capacity to influence its outgoing neighbors (IPO), and the total influence by its incoming neighbors (IPI). The significance of IP stems from the fact that in the evolution of GRNs, particularly in cancer, the mutation of gene functions is usually caused by the localized rewiring of influential genes [Ruan et al., 2015]. To compare two GRNs, the differences in their functions are compared, and a larger functional gap implies larger distance between them. Therefore, paying more attention to high-IP genes will lead to a more effective and credible result. The detailed description of IP is discussed in Section 3.

GHD is a metric that measures the distance of two networks proposed by [Ruan et al., 2015]. Given the adjacency matrices of two networks A and B, the GHD of them is defined as:

where N = |V | is the scale of A and B, a′ij and b′ij are meancentered edge weights which are based on the original edge weights aij and bij.

3 Methodologies

In this Section, the details of FSMmethod is described. The framework of FSM is shown in Figure 2. It can be divided into two parts. First, The influence power from one gene to another is calculated according to their local topological structure. In the second part, each node's IPI and IPO are derived according to a specific rule and assembled into a vector which is termed as total influence power vector (T IV ). Finally, the similarity of the GRNs is outputted. Herein, the calculation of IP and the previously mentioned similarity measurement will be detailed.

3.1 Influence Power

Even though the original GHD can measure the difference between two GRNs, it has certain weaknesses: (1) It is difficult to apply to unweighted GRNs, and thus weights should be manually assigned to each edge. (2) For large-scale GRNs, calculations involving large adjacency matrices should be performed, which requires considerable memory space and CPU resources. Accordingly, IP is proposed to determine efficiently influential genes with the largest contribution to the differences between GRNs. As stated in Section 2.2, the basic principle of IP is to rank genes in a GRN based on their influence during the regulation process. Hence, the most influential genes should be determined. The value of IP of u on v is obtained in terms of the available influence power from u and the influence component from u to v, which are denoted by available power (AP(u)) and influence component (IC(u → v)), respectively. As different neighbors of u have different local topological structure, the intensity of influence is also different. The former measures the "power" of gene u, whereas the latter measures the probability of gene u influencing v.

If u has a large number of incoming neighbors, then u is regulated by several different sources. These incoming neighbors are regarded as the power providers. An influential gene can exert more "power" on other genes. Hereby, the numbers of genes in node sets inN 1 (u) and inN 2 (u) are used to present the influence of the provider.

Definition 1. Given a directed GRNG = (V,E), the available power of a gene u is denoted by AP(u) and is defined

In Equation 2, ∣∣inN 1 (u)∣∣ represents the width of pre-nodes of u, while ∣∣inN 2 (u)∣∣ represents the depth of pre-nodes of u. When considering the available power of a specific gene, its pre-nodes' depth property and width property are equally important.

Considering the gene G2 of GRN A in Figure 2, the volume of "powers" it distributed to its neighbors must not be the same due to the different topological structures of these neighbors. To quantify the intensity with which a gene uses its available power to influence its neighbors, we make an assumption that a gene tends to indirectly influence more other genes in a network. This assumption is very intuitive if we imagine the influence power as a flow in network and flow seeks its own level. For two genes u and v, given two gene sets φ1 = outN 1 (u) + outN 2 (u) and φ2 = v + outN 1 (v), φ1 and φ2 would together determines the chance that gene v would be regulated by u.

Definition 2. Given a directed GRNG = (V,E), the influence component (IC) between two genes u and v is defined as the ratio of their degree distributions, that is,

Finally, by considering these two factors, the IP of gene u on gene v is modeled as follows.

Definition 3. Given a directed GRN G = (V,E), for two genes with regulation relationship (u, v) ∈ E, the influce power (IP) from u to v is defined by

As the differences in GRN functions, the IPs of the same gene across different GRNs are also different. Hence, from the differences of the IPs, the distance between two network can be estimated. Then the FSM method based on the enhanced GHD is introduced to measure this type of distance.

3.2 Implementations of FSMbased on Influence Power

To overcome the second weakness described in Section 3.1, the size of input data should be reduced. A desirable GRN similarity measurement metric should satisfy two requirements: First, it should handle real-world GRNs with acceptable time overhead and spatial cost. Secondly, the distance between two different GRNs should be normally distributed and be able to capture the subtle variations. For the first requirement, a total IP vector (TIV) that represents IPs of all nodes is constructed to replace the original GRN adjacency matrix for dimension reduction. For the second requirement, an enhanced GHD for vectors, called vGHD, is proposed to represent as the similarity.

Definition 4. Given each edge's IP, IPO and IPI are defined as follows:

where p1 ∈ outN 1 (i) and q1 ∈ outN 1 (p1), approximately, q2 ∈ inN 1 (i), p2 ∈ inN 1 (q2).

After we obtain the two vectors, IPO and IPI. TIV is defined as follows:

where S(·) is the softmax function for normalization which smooth the scale gap between IPO and IPI .

Given a GRN, the steps of computing TIV are given in Algorithm 1. TIV is a vector storing each gene's IP. Compared to the raw adjacency matrix, it has smaller size and its elements are the weights of different genes.

Definition 5. Given two GRNs GA, GB and their scale N, their similarity, denoted by vGHD, is defined as follows:

Algorithm 1 TIV Calculation

Input: A GRN G = (V, E)

Parameter: N = |V |

Output: all genes' IPs in G

• 1: Initialize TIV to be a zero vector with 1×N elements.
• 2: Initialize a temp matrix Λ with N ×N elements.
• 3: for gene u in V do
• 4: Let O∗ be the set of outgoing neighbors of gene u.
• 5: for gene v in O∗ do
• 6: Calculate IP(u → v) by Equations (2)-(4)
• 7: Λ(u, v) = IP(u → v)
• 8: end for
• 9: end for
• 10: Calculating IPO based on Λ following Equation (5).
• 11: Calculating IPI based on Λ following Equation (5).
• 12: TIV = S(IPO) +S(IPI)
• 13: return TIV

where

TIVAi is the i-th component of TIV derived fromA. TIVBi is analogously defined.

3.3 Distribution Analysis for vGHD

For measurement purposes, the similarity results should be credible and interpretable. Assuming that there is a typical GRN A from brain cancer, and a new incoming GRN B from an unkonwn person. To determine whether B is also brain cancer (as shown in Figure 1), the distribution of vGHD(A, B) should be analyzed.

First, the null hypothesis is set as follows:

H0 : networks GA and GB are independent.

Under the null hypothesis, the definition of tvGHD in Equation (9) can be rewritten as follows:

where c is a constant. By the transformation of Equation (7), well-known sufficient conditions for asymptotic normality can be used, which can also be easily verified in practice. With the sufficient conditions in [Friedman et al., 1983; Pham et al., 1989], it can be proved that the distribution of vGHD in Equation (9) is approximately normal, that is,

where µ and σ are the mean and standard deviation, respectively. The simulation results shown in Figure 3 demonstrate that the distribution of vGHD is nearly normal distribution.

As vGHD(A,B) obeys follows a normal distribution, given a threshold α, if the p-values exceed α (e.g. 0.05), then H0

cannot be rejected, and thus A and B are independent. Higher p-values indicate a higher probability that B is not brain cancer. Hence, the distribution of vGHD provides strong evidence of FSM's interpretability and credibility.

3.4 Computational Complexity

The proposed method considers the local topology structure of the genes when calculating IP. Its computational complexity depends on the density of the input graph. Let k be the number of average neighbors of a gene, and k2 be the number of second-order neighbors. According to Equations (2)- (5), the computational complexity for calculating each node's IP is O(k2n). As k is constant and significantly less than n, it is in fact nearly O(n). Furthermore, according to Section 2.2, the computational complexity of the original GHD algorithm is O(n2 ). However, after the improvement in the form of Equation 7, it reduces to O(n). The total computational complexity of the proposed method is O(k2n + n) ≈ O(n), k ≪ n.

4 Experiment

4.1 Datasets

Several public datasets of different scales and sources were used; they are all related to gene regulation: AIDS [Riesen and Bunke, 2008], BZR [Sutherland et al., 2003], COX2 [Sutherland et al., 2003], NCI1 [Wale et al., 2008], ENZYMES [Borgwardt et al., 2005], and PROTEINS [Borgwardt et al., 2005]. Furthermore, there are three real-world GRNs derived from gene expression profiles were used, namely, IDH, OC [Zhang et al., 2016], and RAT [Stevenson et al., 2007]. IDH contains 12717 gene regulation relationships from brain tumor data collected from TCGA panglioma samples, and is divided into two subtypes: IDHwildtype and IDH-mutation. OC comprises GRNs of patients with ovarian cancer and consists of 11750 genes. There are two classes in this dataset according to the treatment response (platinum-sensitive and platinum-resistant). RAT comprises GRNs with 8799 genes of two groups of mice: one is the experimental group that was exposed to smoke for about 200 days, and the other is the control group that lived in natural environment.

4.2 Performance Evaluation

In the experiment, the proposed method was compared with these baselines and the following state-of-art methods under the experiment setting described in Section 2.1:

• Randomwalk kernel (RW) [Vishwanathan et al., 2010], it counts common walks in two graphs.
• Shortest path kernel (SP) [Borgwardt and Kriegel, 2005], it counts shortest paths of equal length between pairs of nodes.
• PM [Nikolentzos et al., 2017], it finds an approximate correspondence between the sets of vectors of the two graphs.
• Propagation kernel (PK) [Neumann et al., 2016], a general graph-kernel framework for efficiently measuring the similarity of graphs. It is based on monitoring how information spreads through a set of given graphs.
• Hash graph kernel (HG) [Morris et al., 2016], where continuous attributes are continuously turned into discrete labels using randomized hash functions.

Table 1 shows the average prediction accuracy and standard deviation. Table 2 shows the time required for calculating the distance between two networks.

The experiments were conducted on GRNs and biological networks that, to some extent, have relationships with the process of gene regulatory. The results shown in Table 1 demonstrate two points: (1) In large networks with higher average degree, such as ENZYMES, PROTEINS, and IDH, the proposed method performs the best. The obtained results conforms to the intuition that influential genes should be weighted more in network comparison. Furthermore, the genes' contributions to the mutation of the entire network are different, i.e., influential genes (usually hub genes) contribute more than other marginal genes. (2) In some small networks with lower average degree, although the proposed method is not the best, it still performs well. This is because in those networks, the aggregation effect is not pronounced. Figure 4 shows the comparison of degree distribution between OC, RAT, COX2, and BZR. It can be seen that in OC and RAT, there are only a few genes with large degree distributions that control other nodes. By contrast, in COX2 and BZR, the degree distribution is even across all nodes. The essence of FSM is to detect the differences of those few genes that control the remaining genes. These are assigned higher weights. When the connectivity of higher-position genes changes, a higher distance score is added to them. Hence, the relative accuracy of FSM varies across different datasets.

In terms of runtime, FSM is significantly faster than the other methods. Compared to the RW, FSM is nearly 300 times as fast. This enables FSM to handle large-scale network similarity measurement problems, particularly realworld GRNs such as OC, RAT, and IDH, which usually contain thousands of nodes. In smaller-size networks such as BZR and COX2, FSM also performs well in terms of computation time. However, in AIDS and PROTEINS, PK is a little faster than FSM. It seems that the speeds of PK and FSM are tied in small datasets, but FSM is more suitable in handling larger data just like IDH, OC and RAT.

5 Conclusion

The FSM method was proposed for fast similarity measurement between GRNs via vGHD and influence power. A new index (IP) was defined to describe influential genes in the network. Inspired by common phenomena, the basic principle of FSM is that more influential genes contribute more to the differences of the GRNs. The proposed method exhibits good performance in GRNs and other related bio-networks.

Acknowledgements

This work was supported in part by the National Natural Science Foundation of China (No. 61572369), in part by the Key Projects of Guangdong Natural Science Foundation (No. 2018B030311003)