An approach based on Boolean dynamics of biological networks for disease candidate gene prioritization  

Duc-Hau Le

Abstract

One of the challenging issues in biomedicine is to identify candidate genes likely to be associated with specific diseases. To tackle this issue, many computational methods have been developed for ranking disease candidate genes according to their relevance to a disease of interest. This relevance has been defined as functional similarity or relatedness between candidate genes and known disease-associated genes. In addition, it is widely accepted that most functions of cell are controlled by biological networks constituted by interactions between genes, proteins, metabolites and other molecules. Defects by mutations on genes/proteins may cause a disease occurred in a person. Also, these mutations may affect other genes/proteins through structures of the biological networks. In this study, we propose a novel method, namely BoolDGP, to assess the relevance of candidate genes to a disease of interest by measuring the degree of mutational effect from known disease-associated genes to candidate genes. Particularly, we mutate known disease-associated genes and measure the effect of this mutation on candidate genes based on Boolean dynamics of biological networks. Based on this value, candidate genes can be prioritized and finally top-ranked candidate genes can be selected as novel promising disease genes. Simulation results on a set of diseases showed that the proposed method is superior to a state-of-the-art one, which is based on a random walk with restart algorithm. Using the proposed method, we identified 27 genes associated with breast cancer with evidences from literature.

Introduction

Disease gene prioritization, the task of predicting the most plausible candidate disease genes, is an important issue in biomedical research, and a variety of approaches have been proposed (1, 2). Identification of disease-associated genes also leads to more effectively researches about therapies for genetic diseases (e.g., hypertension, breast cancer, etc…) and gradually approaches a future of personalized medicine. In past decades, linkage analysis was usually used to identify novel disease genes, in which susceptible loci including hundreds of genes are investigated, and thus it is much costly for doing many experiments in wet lab. Therefore, ranking/prioritization methods for such candidate genes are introduced (genes are ordered by their relevance to a disease of interest). Highly ranked genes are further investigated to find out associated biomedical evidences. Therefore, the goal of gene ranking/prioritization is to identify novel disease-associated genes.

The candidate gene prioritization methods have been proposed by three main directions: i) Functional annotation-based; ii) Network-based; and iii) Machine learning-based. In which, functional annotation-based methods have prioritized candidate genes by measuring the degree of similarity of each candidate genes to a set of known  disease genes based on profiles which were built from many functional  annotation data sources (3, 4, 5). Therefore, those methods mostly focused on the integration of various biological datasets in order to obtain more accurate similarity. However, those approaches are limited in that functional annotation data sources have not covered whole human genome yet. Beside the functional annotation-based methods, machine learning-based techniques based on binary classifiers for identification of disease-associated genes have been also studied. At the early, machine learning-based studies usually approached disease gene prediction as a binary classification problem. Indeed, a number of binary classification techniques have been proposed to the problem such as decision tree (6, 7) k-nearest neighbors (8), Naïve Baysian classifier (9, 10), artificial neural network (11, 12) and support vector machine (13, 14, 15). In these studies, the learning samples are comprised of positive training samples and negative training samples. In which, positive training samples are constructed from known disease genes, whereas negative training samples are the remaining which are not known to be associated with diseases. This is the limitation of binary classifier-based solutions for disease gene prediction problem, since the negative training set should be actual non-disease genes. However, construction of this set is nearly impossible in biomedical researches because it is often the case in biomedicine that not observing an association does not imply the association does not exist (there are no proven negatives). In fact, the unknown set may contain unknown disease genes; therefore, to reduce this uncertainty of the earlier methods, semi-supervised methods have been proposed to the problem (16, 17, 18), where the classifier is learned from both labeled (i.e., known disease genes) and unlabeled (i.e., the unknown genes) set. However, negative samples must be still defined in those studies. To overcome limitations of these two main approaches, network-based methods for identification of disease-associated genes have been proposed (19, 20). These network-based methods mainly based on biological networks and therefore they are not limited by the coverage of functional annotation data sources. In addition, they dominate the other two main approaches is because they are based on a principle of “disease module” (e.g., genes/proteins associated with the same or similar diseases usually form functional/physical modules on gene/protein interaction networks (21, 22, 23)). Furthermore, network-based methods target to the essence of disease gene prediction problem, which is ranking/prioritization, instead of classification of candidate genes (i.e., assign a candidate gene to either disease genes or non-disease gene class) as machine learning-based methods do.

It is widely accepted that defects by mutations on genes/proteins may cause a disease occurred in a person. Also, these mutations may affect other genes/proteins through structures of the biological networks. This phenomenon has an underlying reason that mutations on genes/proteins affect the robustness of biological networks and propagation of the effect is controlled by structural properties of biological networks. In this study, we propose a new method, namely BoolDGP, based on Boolean dynamics of biological network to measure the effects from known disease genes to candidate genes. Based on this value, candidate genes can be prioritized and finally top-ranked candidate genes can be selected as novel promising disease genes for further experimental studies. Comparing performance of the proposed method on a set of 25 disease phenotypes with that of a state-of-the-art network-based method, Random Walk with Restart (RWR) algorithm (24), we found that our method outperforms the RWR-based one. In addition, we used the proposed method to identify novel breast cancer-associated genes. Interestingly, 27 out of 50 highly ranked candidate genes are evidenced to be associated with breast cancer.

Materials and Methods

Biological networks and known disease gene associations datasets

To assess mutational effects of known disease-associated genes to candidate genes, we used a large-scale human signaling network collected from a published study (25). This network consists of 1,539 nodes and 4,754 interactions and is the largest one available on literature. In addition, a set of diseases and their known associated genes were collected from OMIM (26). Because only small set of genes are available on the current human signaling network and only diseases with at least two known associated genes are suitable for leave-one-out cross validation analysis; therefore, a set of 25 disease phenotypes was finally selected.

Boolean dynamics-based measure for prioritizing candidate genes

To define a measure assessing the mutational effects of known disease-associated genes on candidate genes, we employed a Boolean network model, which has been widely used to represent biological networks and successfully captured some biological characteristics (27, 28, 29, 30, 31, 32, 33). In particular, it has been also frequently used in simulating the dynamics of various signaling networks (34, 35, 36, 37, 38, 39, 40, 41, 42).

A random Boolean network

When a Boolean network is represented by a directed graph G (V, E), each v_i∈V  has a value of 1 ("on") or 0 ("off"), which represents the possible states of the corresponding elements. The value of each variable v_i at time t+1 is determined by the values of k_i other variables with a link to v_i at time t by the Boolean function. Hence, we can write the update rule as   where we randomly select either a logical conjunction or disjunction for all signed relationships in f_i with a uniform probability distribution. For example, if a Boolean variable v has a positive relationship from v₁, a negative relationship from v₂, and a positive relationship from v₃, then the conjunction and disjunction update rules are (t) and (t), respectively. In the case of a conjunction, the value of v at time t + 1 is 1 only if the values of v₁, v₂, and v₃ at time t are 1, 0, and 1, respectively whereas, in the case of a disjunction, the value of v at time t + 1 is 1 if at least one of the states of the clauses, v₁(t), , and v₃(t) is 1. Although there can be many other logical functions in addition to conjunction and disjunction functions, biological networks were successfully described by Boolean models using only those two functions in many previous studies (31, 32, 33, 43, 44, 45, 46). In addition, the sign of each link is determined between positive and negative ones uniformly at random.

Given a Boolean network with N Boolean variables, v₁, v₂, …, v_N, we define a network state as a vector consisting of values of the Boolean variables: there are 2^N states in total. Each state transits to another state through a set of N Boolean update functions, f₁, f₂, …, f_N. We can construct a state transition diagram that represents the transition of each state. A state trajectory starts from an initial state and eventually converges to either a fixed-point or a limit-cycle attractor. Attractors can represent diverse behaviors of biological networks, such as multi-stability, homeostasis, and oscillation (47, 48, 49). In addition, we define a transient sequence of values of a node v as follows: When a Boolean network G(V, A) was initialized with v₁(0), v₂(0), …, and v_N(0) at the starting time 0, v_i(t_0,t₁) represents a sequence of the transient values of a node v_i during the time interval from t₀ to t₁.

Effectiveness from a node to another node in arandom Boolean network

In Boolean networks, we propose a novel measure, effectiveness, to quantify the influence from a node to anothernode in terms of the network dynamics. To define it, we first introduce two types of perturbations, an initial-state perturbation and an updating-rule perturbation. Given a Boolean network initialized with v₁(0), v₂(0), …, and v_N(0), the initial-state perturbation at a node v_iÎV means flipping v_i(0) to . On the other hand, the updating-rule perturbation at a node v_iÎV means switching the updating-rule at v_i from a conjunction function to a disjunction function or vice versa, depending on the current function type. Assuming a perturbation at v_i, we define the effectiveness from v_i to another node v_j, e(v_i, v_j), as follows:

1. Let τ_i, the valid convergent time of v_i, defined as τ_i = max{T_i, T’_i} where T_i or T’_i represent the time steps for the network to converge to an attractor when v_i was subject to the perturbation or not, respectively.
2. We obtain two different transient sequences of v_j, v_j(0,τ_i) and v’_j(0,τ_i), when v_i was subject to the perturbation or not, respectively.
3. Then, we compute e(v_i, v_j)=d(v_j(0,τ_i), v’_j(0,τ_i)) / τ_i where d(∙) means the Hamming distance (i.e., the number of bits having different values) between two sequences. Thus, e(v_i, v_j) represents how largely the trajectory with respect to v_j was affected by the perturbation at v_i. This also measures the mutational effect of v_i to v_j.

Fig. 1. An illustrative example of calculating effectiveness in a Boolean network. A Boolean network with 8 nodes and 14 links where arrows and bar-headed lines represent positive and negative interactions, respectively. ‘AND’ and ‘OR’ denote conjunction and disjunction update functions, respectively. Trajectories starting from an initial state (11010010) and another state (11011010) where v₄ is subject to an initial-state perturbation. States of the network (i.e., eight-bit strings in rectangles) represent values of v₀ through v₇ in sequence and grayed rectangles with dashed lines mean attractors. These states and trajectories are calculated from the network. Then, effectiveness from v₄ to v₇ with respect to the trajectories was computed

Therefore, effectiveness is a measure about how largely each node is affected by perturbation at the other node in terms of dynamics. In a Boolean network, a node is called a functional important node if a perturbation at the node makes the network converge to another attractor, which is different from the original attractor to which the network converged when the node was not subject to the perturbation. In this regard, disease genes can be considered as important nodes in signaling networks and the effectiveness in Boolean networks can be used to represent the effectiveness on candidate genes when known disease-associated genes are mutated in signaling networks. Figure 1 shows an example of the calculation of effectiveness of v₄ to v₇. To compute e(v₄, v₇), we get two transient sequences of v₇, v₇(0,τ₄) and v₇’(0,τ₄), when v₄ was subject to a perturbation or not, respectively.

Effectiveness from a set of nodes to a node in arandom Boolean network

In a similar way, given an initial state, the effectiveness from a set of nodes S to a node can be calculated by applying perturbations on all nodes in the set simultaneously, and then the effectiveness from these nodes to a node v_je(S, v_j) is calculated. For a set of initial state I_s, the effectiveness from a set of nodes to a node is formally defined as following:

For the identification of disease-associated genes, S is a set of known genes of a disease of interest, v_i is one of candidate genes, then  measures the effectiveness from known disease genes to a candidate gene. Therefore, candidate genes can be ranked to be associated with the disease of interest by this measure.

A. Random Walk with Restart (RWR) algorithm

To show the advance of our proposed method, we selected a state-of-the-art network-based method, which is based on a random walk with restart (RWR) algorithm. RWR is a variant of the random walk (50) and it mimics a walker that moves from a current node to a randomly selected adjacent node or goes back to source nodes with a back-probability gÎ(0, 1). RWR can be formally described as follows:

where P^t is a N´1 probability vector of |V| nodes at a time step t of which the ith element represents the probability of the walker being at node v_i∈V, and P⁰ is the N´1 initial probability vector where the value of an element corresponding to a non-source node or a source node is zero or 1/|S|, respectively. S is the set of source nodes. The matrix W’ is represented by a transition probability matrix and thus (W’)_ij, the (i, j) element in W’, denotes a probability with which a walker at v_i moves to v_j among V\{v_i}. Formally, for an unweighted network, it is defined as follows:

where (V_out)_i is a set of outgoing nodes of v­_i.

Fig. 2.  Comparison of performance between BoolDGP and RWR. The performance of two methods was assessed using LOOCV method on the set of 25 disease phenotypes from OMIM. For BoolDGP, initial-state perturbation method and a set of 100 initial states were used. For RWR, the back-probability was set to 0.5.

All nodes in the network are eventually ranked according to the steady-state probability vector P^¥, which is obtained by repeating the iterations until ||P^t+1-P^t||<10^-6 in this study.

For the identification of disease-associated genes, S is a set of known genes of a disease of interest, an element of P^¥measures how much relevant to S a gene in the network is. Therefore, it is the degree of association between a candidate gene and the disease of interest.

Performance Evaluation

Ranking performance was assessed through the leave-one-out cross-validation (Shortly called LOOCV) process. For each disease phenotype (d), in each round of LOOCV, we held out one known d-associated gene. The rest of known d-associated genes are specified to a set of source nodes (i.e., S). The held-out gene and remaining genes in the human signaling network, which were not known to be associated with d, were ranked by the two ranking methods. Then, we plotted the receiver operating characteristic (ROC) curve and calculated the area under the curve (AUC) to compare the performance of these two methods. This curve represents the relationship between sensitivity and (1-specificity), where sensitivity refers to the percentage of known disease-associated genes that were ranked above a particular threshold and specificity refers to the percentage of genes which were not known to be associated with a disease ranked below this threshold.

Results

Performance Evaluation

To assess performance of BoolDGP in ranking disease candidate genes, we randomly selected I_s = 100 initial states and used LOOCV to draw the ROC curve and calculate AUC value for the set of 25 disease phenotypes collected from OMIM (See Materials and Methods section). To compare the performance of our method with that of RWR-based method, we follow the same procedure as did for BoolDGP, and set the back-probability to 0.5 since the performance of the RWR-based method is stable with the changing of the back-probability parameter (24). Figure 2 shows the performance of the two methods, it is clear that BoolDGP (AUC = 0.77) outperformed the RWR-based method (AUC = 0.73). This result implied that the disease gene prioritization problem can be approached effectively with the hypothesis that the mutation on known disease genes of a disease of interest affect other genes and this effectiveness can be used as the degree of association between candidate genes and the disease of interest.

Case study: Breast cancer

In this section, we show the ability in predicting novel genes associated with specific disease. In particular, we test our method with breast cancer (OMIM ID: 114480), this is a complex disease with 22 known associated genes. However, only ten of them are available on the human signaling network. Using these known genes as source nodes and considering other genes in the network as candidates, we calculated the effectiveness from these known genes to all the candidate genes. Candidate genes then are ranked based on the effectiveness. We selected 50 highly ranked candidates to find evidences of their association with breast cancer from literature. Interestingly, a total of 27 among them have at least one evidence to be associated with breast cancer (Table I). For example, a study showed that genomic rearrangements of the CHK1 gene is associated with breast cancer (51). Also, polymorphisms/haplotypes in GADD45A contribute to breast cancer risk, at least to sporadic breast cancer (52). In addition, high expression of HBEGF is related to the biological aggressiveness of breast cancer (53). mTOR is a critical target for survival signals generated by phospholipase D in breast cancer cells (54). Other candidate genes in the top 50 can be promising ones for future studies.

Table 1: A total of 27 among top 50 ranked candidate genes are evidenced to be associated with breast cancer from literature

Conclusions

In this study, we proposed a novel network-based method for identification of disease-associated genes. This is based on Boolean dynamics of biological networks under a hypothesis that mutation on known disease genes of a disease of interest affect other genes through the network and this effectiveness can be used as the degree of association between candidate genes and the disease of interest. Simulation results showed that our method is superior to a state-of-the-art network-based method. Using the proposed method, we also predicted 27 novel breast cancer-associated genes with evidences from literature. In future study, besides the Boolean dynamics, we are going to integrate structural properties of biological networks to a measure of association between candidate genes and diseases of interest since some studies have shown special structural properties of disease genes in the biological networks.

Acknowledgement

This work was supported by Vietnam Institute of Advanced Study of Mathematic (VIASM), Ministry of Education and Training under contract number 76NC/2014/VNCCCT.

Reference