Problem formulation
Let H be a bipartite graph consisting of two disjoint sets of vertices, where one set represents all known human genes {g1, g2, . . . , g
N
}, while the other set represents all known genetic diseases {d1, d2, . . . , d
r
}. The associations between those genes and genetic diseases can be obtained from either the Online Mendelian Inheritance in Man (OMIM) database [22] or similar databases.
Although a disease d
k
may associate with one or several genes, the number of all known disease genes m is much smaller than N . Hence, associations of most other genes are still not known and need to be analyzed. Without loss of generality, we can reorder the set of all human genes as a vector (g1, g2, . . . , gN ), according to a set of given gene-disease associations, where gn+1, gn+2, . . . , gn+mare genes associated with at lease one known disease (disease genes), and g1, g2, . . . , g
n
are others. Here N = n + m, and n is the number of all genes that are not known to associate with any diseases and they are called unknown genes in this paper.
For a specific disease d
k
, the issue of disease gene identification is to find a set of candidate genes which may have associations with d
k
. To achieve this, let be a vector of binary class labels (i.e. taking the value zero or one) defined on all genes, where represents gene g
i
being a disease gene of d
k
, and otherwise. Since we have to address each genetic disease one by one, we take d
k
for example, and ignore the superscript k in the vector xkfor simplicity as x = (x1, x2, . . . , xn+m), hereafter. Therefore, the identification of disease genes is equivalent to find labels of x
i
for all unknown genes. Identification of disease genes for other diseases can be similarly conducted by changing d
k
to another disease.
In this paper, the issue of disease gene identification is formulated as a two-class classification problem by using Bayesian analysis and logistic regression. The conditional probability p(x
i
= 1| Φ) for each unknown gene is first calculated in an inference stage, and a decision score is then obtained according to this probability in a decision stage [23]. Here Φ represents the information used to make the inference, such as a vector of prior labels of x, the connectivity of the bipartite graph H, the neighborhood relationships of g1, g2, . . . , g
N
, and the similarity relationship between d1, d2, . . . , d
r
. The flow diagram of the proposed algorithm is depicted in Figure 1.
Prior label estimation
The logistic regression based algorithm needs a vector of prior labels for x. For those known disease genes, one can directly assign 1 or 0 according to the known gene-disease associations. For those unknown disease genes, a prior probability of each gene get the label 1 should be first estimated.
The simplest way is to assign the prior probability as 0 for all unknown genes. The prediction results in this case is denoted as P0 hereafter.
However, one can make it better by using additional prior information, such as protein complex data, to estimate prior probabilities for unknown genes. This is not only due to the fact that they are naturally available from various databases, but also due to their capability to describe the module characteristic of disease genes. If a disease is resulted from the disfunction of a protein complex, then any component of the complex should associate with the disease with a high probability.
Similar to the method used in [19, 20], if a gene g
i
encodes a protein in a complex, then let
be its prior probability, where A is the number of disease genes of the specific disease in the complex, and B is the number of all disease genes in the complex. If g
i
appears in multiple protein complexes, we use the maximum value as its prior probability. If g
i
does not belong to any protein complex, let
be its prior probability, where C is the number of all currently known disease genes of the specific disease, and D is the total number of genes in human genome.
Once a prior probability is estimated for g
i
, its prior label of can be obtained as follows. First, generate a random number following the standard uniform distribution. If the value of the random number is large than , then assign 0 as the prior label for g
i
. Otherwise, assign 1 as the prior label for g
i
. Repeat this step for all unknown genes, one can obtain prior labels for all of them. The prediction results generated by using those prior labels is denoted as P
c
hereafter.
Logistic regression
For a two-class classification problem, each gene is labelled with either 1 or 0. A vector of all binary values of x is called a configuration. In the previous MRF algorithm [19, 20], the configuration x is formulated as an MRF which follows a Gibbs distribution. However, the Markovianity characteristic of the MRF model makes it only considering direct neighbors to construct feature vectors, which limits the capability of the method to use other topological attributes in a biological network. It is also very time consuming to maintain Markov chains for all unknown genes. In [21], we have introduced a logistic regression based algorithm to directly estimate the configuration by using the same feature vectors. However, the application of the logistic regression based algorithm is still limited to single biological network. A multiple data integration method should be further investigated. To generalize the formulation of feature vectors by using other topological attributes and extend its applicability to multiple data integration, we propose an improved logistic regression based algorithm in this study as follows.
Let C1 be a set of genes with label 1 and C0 be a set of genes with label 0. Suppose the following four kinds of probabilities are given: the class-conditional densities p(x|C1) and p(x|C0), which indicate the probability of the configuration x conditional on C1 and C0, respectively, and the class prior densities p(C1) and p(C0), which indicate the prior probability of genes in C1 and C0 being labelled with 1 and 0, respectively.
According to the Bayes' rule, the posterior probabilities of those genes in C1 that are labelled with 1 can be described as a logistic sigmoid function [23, 24]
(3)
and the posterior probabilities of those genes in C0 that are labelled with 0 can be similarly written as
(4)
where the variable t is defined as
(5)
which is related to the four kinds of probabilities.
Although t is often unavailable for a real problem, under general assumptions [23], t can be formulated as a function of different features t = f (·) associated with the integrated networks. To be more specific, let be a prior configuration of all human genes and f be a function. For any given gene g
i
, let ϕ
i
be the feature vector of g
i
that is related to the prior configuration . The posterior probability that the specific gene g
i
has label 1 and 0 are
(6)
and
(7)
respectively. Note that the sum of these two probabilities (6) and (7) must equal to 1 in this two-class classification problem. A linear function f (ϕ
i
) = wT ϕ
i
with variables (feature vectors) ϕ
i
and coefficients (parameters) w is the most commonly used function to ensure the calculation of the posterior probability not too complex.
The key step of the proposed algorithm is the construction of feature vectors. In the previous methods [19, 21], the numbers of direct neighbors that connects to disease genes and non-disease genes are employed as the feature vector for each gene. Take g
i
for example, its feature vector can be written as
(8)
where ϕi1 and ϕi0 are the number of direct neighbors of g
i
that connected to vertices with labels 1 and 0, respectively. It is a three dimensional vector, where the first element represents the constant term. All feature vectors of individual genes together form a feature matrix as
(9)
where N is the number of all human genes. The corresponding parameters are w = (w
0
, w
1
, w
2
)T. Predictions generated by using (9) are denoted as F1 hereafter.
In this study, two extended feature vector construction methods are proposed as follows. Firstly, in a single biological network, not only the number of direct neighbors of g
i
, but also the number of its second order neighbors are employed to construct the feature vector as
(10)
where ϕi1 and ϕi0 are the numbers of direct neighbors of g
i
connected to vertices with labels 1 and 0, respectively, and and are the numbers of the second order neighbors of g
i
connected to vertices with labels 1 and 0, respectively. The contribution of those indirect neighbors has been investigated for predicting disease genes in [20, 25, 26]. The feature matrix in this situation can be written as
(11)
The corresponding parameter vector w = (w0, w1, w2,w3, w4)T is a five dimensional vector. Predictions generated by using (11) are denoted as F2 hereafter.
Secondly, in the multiple data integration situation, suppose there are l biological networks. Let be the number of direct neighbors of g
i
connected to vertices with labels 1 and 0 in the jth network, respectively. The feature vector obtained from those l networks
(12)
is a 2l + 1 dimensional vector. All those feature vectors together form a feature matrix as
(13)
The corresponding parameter vector w = (w0, w1, w2,. . . , w2l−1, w2l)T is a 2l + 1 dimensional vector, and N is the number of all human genes. Predictions generated from (13) by integrating multiple networks is denoted as F3 hereafter.
Parameter estimation
Parameter estimation can be conducted on a training set consists of known disease genes, where known genes associated with d
k
are labelled with 1 and known genes associated with other diseases are labelled with 0. However, as we discussed in [19, 21], the exclusion of most unknown genes reduces the number of vertices with label 0 significantly, thereby making the estimation of parameters inaccurate. Predictions from those inaccurate parameters are unreliable in disease gene identification.
It is noteworthy that the majority of human genes should not be disease genes associated with d
k
. Hence, the inclusion of all unknown genes with prior labels as the training set will make the training set more reasonable, where the number of vertices with label 0 is significantly increased, while the number of vertices with label 1 does not change too much. Such a training set, which consists of both known genes and unknown genes, has proved its powerful and efficient to estimate meaningful parameters in [19–21].
Given a prior configuration for all vertices, a maximum-likelihood estimation (MLE) method can be employed to estimate the parameter vector w. The likelihood function can be written as
(14)
where x
i
is the label of g
i
, ϕ
i
is its feature vector that is calculated according to , f is a linear function of ϕ
i
with the form as f (ϕ
i
) = wT ϕ
i
, and N is the number of all human genes. The log likelihood of (14) is
ln L (w; x1, x2, . . . , x
N
)
(15)
The log likelihood (15) is a convex function [27]. Hence, we can find an unique global optimal solution by solving a convex optimization problem. In this study, the standard MATLAB function fminunc() is employed to find a numerical solution of (15) (by calculating the minimum of − ln L (w; x1, x2, . . . , x
N
)). The the initial value of w is simply set as zero for the fminunc() function.
Decision score and evaluation methods
The logistic regression based algorithm returns a set of posterior probabilities during the inference stage. One can directly use those probabilities to make decisions in the following decision stage. However, the posterior probabilities do not always work well due to the hardness to set a threshold for a genetic disease. Inspired by the DIR method [18], we propose to use a percentage value of a posterior probability as the decision score for each gene. The decision score is calculated as follows
(16)
where {p1, p2, . . . , p
n
} is the posterior probabilities of individual unknown genes, and q
i
is the top percentage value of p
i
among all those posterior probabilities. A candidate gene is more likely to be associated with d
k
, if its decision score is larger than majority of others.
To evaluate the performance of the proposed algorithm, the leave-one-out cross validation paradigm is employed by using above decision scores. The receiver operating characteristic (ROC) curve is employed as one of the evaluation criteria, which shows the relationship between the true positive rate (TPR) and the false positive rate (FPR) by varying a threshold for determining positives. The area under the ROC curve (AUC) is employed to show the overall performance of algorithms.
The positive control genes are those known disease genes associated with d
k
. For those negative control genes, although they are indispensable to calculate false positives and true negatives, it is generally hard to obtain a true negative dataset [28]. In this study, the negative control genes are randomly selected from known disease genes that do not associate with d
k
. Since those genes have been widely studied as disease genes for other genetic diseases, it is less likely for them to be disease genes for a different specific disease. If there are s known disease genes associated with d
k
, we randomly select such genes as a negative control set. Each gene belonging to the negative control set is also validated by using the leave-one-out cross validation paradigm.
The proposed algorithm is compared with four previous algorithms: (1) the initial logistic regression based algorithm proposed in [21]; (2) the RWR algorithm proposed in [17]; (3) the MRF algorithm proposed in [19]; and (4) the DIR algorithm proposed in [18]. The first algorithm is applicable to a single network. The second and the third algorithms are applicable to both single network and multiple data integration. The fourth algorithm works only for multiple data integration. All those algorithms identify disease genes with high prediction performance and they work better than many previous methods [17–19, 21].
Algorithm
The step-by-step description of the proposed logistic regression based algorithm is given as follows.
Input: The vector of all human genes (g1, . . . , gn+m), where (g1, . . . , g
n
) are unknown genes, and (gn+1, . . . , gn+m) are known genes; l integrated biological networks G1, G2, . . . , G
l
; a set of protein complexes; and a set of gene-disease associations.
Output: The vector of decision score for each unknown gene for each disease.
-
1:
For a specific disease d
k
, calculate prior probabilities for all human genes, where the prior probability of unknown genes are calculated according to (1) and (2).
-
2:
For each known gene gn+i, i = 1, . . . , m, if gn+iis known to be associated with d
k
, let Otherwise, let
-
3:
Assign prior labels for all genes according to the prior probabilities respectively.
-
4:
Calculate the feature vector ϕ
i
for each g
i
according to the integrated biological networks and .
-
5:
Estimate parameters by maximizing the log likelihood ln in (15) based on and ϕ
i
, i = 1, . . . , n + m. A binary logistic regression is performed here by taking the vector as the categorical dependent variables and those label-related feature vectors ϕ
i
as predictor variables. Here i = 1, . . . , N.
-
6:
Calculate the posterior probability p1, . . . , p
n
for each unknown gene according to (6) by using and ϕ
i
.
-
7:
Calculate the decision scores q1, . . . , q
n
according to (16).
-
8:
Repeat all the steps for another disease until every disease is checked.