Graph regularized non-negative matrix factorization with prior knowledge consistency constraint for drug–target interactions prediction

Abstract

Background

Identifying drug–target interactions (DTIs) plays a key role in drug development. Traditional wet experiments to identify DTIs are expensive and time consuming. Effective computational methods to predict DTIs are useful to narrow the searching scope of potential drugs and speed up the process of drug discovery. There are a variety of non-negativity matrix factorization based methods to predict DTIs, but the convergence of the algorithms used in the matrix factorization are often overlooked and the results can be further improved.

Results

In order to predict DTIs more accurately and quickly, we propose an alternating direction algorithm to solve graph regularized non-negative matrix factorization with prior knowledge consistency constraint (ADA-GRMFC). Based on known DTIs, drug chemical structures and target sequences, ADA-GRMFC at first constructs a DTI matrix, a drug similarity matrix and a target similarity matrix. Then DTI prediction is modeled as the non-negative factorization of the DTI matrix with graph dual regularization terms and a prior knowledge consistency constraint. The graph dual regularization terms are used to integrate the information from the drug similarity matrix and the target similarity matrix, and the prior knowledge consistency constraint is used to ensure the matrix decomposition result should be consistent with the prior knowledge of known DTIs. Finally, an alternating direction algorithm is used to solve the matrix factorization. Furthermore, we prove that the algorithm can converge to a stationary point. Extensive experimental results of 10-fold cross-validation show that ADA-GRMFC has better performance than other state-of-the-art methods. In the case study, ADA-GRMFC is also used to predict the targets interacting with the drug olanzapine, and all of the 10 highest-scoring targets have been accurately predicted. In predicting drug interactions with target estrogen receptors alpha, 17 of the 20 highest-scoring drugs have been validated.

Background

According to the DrugBank database, there are over 500,000 drugs for different diseases. However, there are still many diseases for which we have no effective drugs, and there is a long way for drug discovery and drug repurposing. The process of drug discovery and drug repurposing takes some important steps such as finding valid target proteins and discovering proper chemical compounds to interact with the targets (i.e. identifying drug–target interactions or binding affinity between drugs and proteins) [1]. Determining drug–target interactions (DTIs) via wet experiments is both time-consuming and expensive [2, 3]. To increase the probabilities of discovering new drugs or new applications of approved drugs, accurate computational methods to predict DTI are in urgent need to choose a small number of compounds for the wet experiments.

Based on the crystal structure of the target binding site, Cheng et al. [4] constructed maximal affinity model using docking simulations and calculated the maximal affinity values associated with the drug. Campillos et al. [5] used side-effect similarity of drugs to infer the probability of two drug sharing a target. The above two methods could not predict DTIs for targets without known crystal structure or drugs without known side-effects.

To effectively predict DTIs in large scale, a lot of computation models have been introduced. For example, Yamanishi et al. [6] computed drug chemical structure similarities, amino acid sequence similarities of target proteins, and proposed a bipartite graph learning method to predict DTIs by eigenvalue decomposition based on and known drug–target interactions. Based on the same drug structure similarities, sequence similarities of proteins and known drug–target interactions, Bleakley and Yamanishi [7] proposed bipartite local models (BLM) to predict target proteins for a given drug, and the drugs targeting a given protein. Laarhoven et al. [8] proposed Gaussian interaction profile (GIP) kernel method to calculate the similarities between targets (and drugs), and regularized least squares (RLS) was used to predict unknown DTIs. The Kronecker RLS model [9] was proposed based on RLS. However, these methods can not predict DTIs for the drugs or targets if they have no known interactions. In order to compensate the lack of interaction information, based on GIP, Laarhoven and Marchiori [10] proposed a weighted nearest neighbor method to predict DTIs. Mei et al. [11] derived training data from neighbors of new drug (target) candidates, and proposed BLM-NII method by integrating neighbor-based interaction-profile inferring into BLM for DTI predictions.

Instead of utilizing the attributes of drug chemical structure and sequence of target proteins separately, more and more methods integrated multiple features and then built classifiers to make prediction. Wang and Zeng [12] transformed the DTI prediction problem into a two-layer graphical model, used restricted Boltzmann machine to predict diverse types of DTIs, such as direct and indirect interactions. Based on multiple information, random walk with restarts (RWR) was used for DTIs feature extraction on positive-unlabeled learning method [13] and DTINet method [14]. In order to predict large scale DTIs, convolutional neural network (CNN) was applied to to extract drug and target features information [15,16,17,18]. FCS [19] and KGE_NFM [20] were proposed to obtain the low-rank representations for information of multi-omics in DTI prediction. Other information fusion methods including negative sample screening framework [21], SITAR framework [22] and multiple kernels learning [23, 24] were successively proposed. However, the computation complexity of large scale similarity matrices is high.

The known DTIs are usually represented by a drug–target interaction matrix, and matrix factorization method has been widely used in DTI prediction. Matrix factorization is a embedding model that is used to decompose interaction matrix into two feature matrices of low ranks which represents the interactions between drugs and targets as the product of feature matrices. For example, Gönen [25] proposed a kernelized Bayesian matrix factorization with twin kernels method to predict DTIs. The twin kernel matrices were constructed by chemical similarity function of drug compounds and genomic similarity function of target proteins. Combined logistic matrix decomposition with neighborhood regularization [26] and a variational Bayesian multiple kernel logistic matrix factorization method [27] were proposed to infer interactions. MSCMF [28], GRMF [29], $L_{2,1}$-GRMF [30], GRGMF [31], HCNMF [32], CHNMF [33] and SRCMF [34] were proposed by means of graph regularization. However, their performances are still not satisfying due to the lost information of these decompose strategies.

To implement matrix factorization, the above methods used either the alternating least squares algorithm [35] or the multiplicative update algorithm [36]. However, there are other efficient algorithms to implement matrix factorization. Hoyer [37] introduced a sparseness measure, and proposed a projected gradient descent algorithm for non-negative matrix factorization with sparseness constraints. Lin [38] proposed an improved projected gradient method for NMF with bound constrains. It has not been proved that the alternating least squares algorithm and the multiplicative update algorithm converge to a stationary point when convergence does occur [39]. Due to the use of a simple geometric rule for the step size, gradient descent methods are very sensitive to the initialization [40] and often produce a poor factorization. To overcome the limitations of the above methods, the alternating direction method (ADM) has attracted attention. To obtain a higher-quality matrix factorization with less computing time, Zhang [39] extended the alternating direction method (ADM) to solve NMF. Based on the above algorithm, Xu et al. [41] devised an improved alternating direction algorithm (ADA) to solve the non-negative matrix factorization-and-completion problem.

The high-dimensional data are in fact sampled from a nonlinear low-dimensional manifold embedded in the high-dimensional space, and according to [42], the model learning performance can be greatly improved if the intrinsic geometrical structure of the manifold have been taken into account. Shang et al. [43] showed that a graph dual regularization helps NMF improve clustering performance since the graph dual regularization considers the underlying geometric structures of both the data manifold and the feature manifold. In this paper, in order to predict DTIs more accurately and quickly, we propose an Alternating Direction Algorithm to solve Graph Regularized non-negative Matrix Factorization with prior knowledge consistency constraint (ADA-GRMFC). The prior knowledge consistency constraint aims to ensure that the decomposition result is consistent with the prior knowledge of known DTIs. The alternating direction algorithm ensures that ADA-GRMFC can converge to KKT point.

Extensive experimental results show that ADA-GRMFC has better performance than other state-of-the-art methods. In case studies involving the drug olanzapine and the target estrogen receptor alpha, all the 10 highest-scoring targets predicted to interact with olanzapine, and 17 of the 20 highest-scoring drugs predicted to interact with estrogen receptor alpha have been validated by wet experiments. The case studies show that, for drugs that do not have any known target proteins and for proteins that are so far not approved as drug targets, ADA-GRMFC also has good prediction performance.

Materials

The experimental data include known drug–target interactions, drug chemical structures and target protein sequences. They were from public databases BRENDA [44], KEGG BRITE [45], SuperTarget [46] and DrugBank [47] and were downloaded from the website: http://web.kuicr.kyoto-u.ac.jp/supp/yoshi/drugtarget/.

The target proteins have the following four types: nuclear receptors (NR), G protein-coupled receptors (GPCR), ion channels (IC) and enzymes (E). Accordingly, the benchmark known drug–target interactions are divided into four datasets NR, GPCR, IC and E. The sizes of the four datasets are different. In the NR dataset, there are 90 known interactions between 54 drugs and 26 nuclear receptors; in the GPCR dataset, there are 635 known interactions between 223 drugs and 95 G protein-coupled receptors; in the IC dataset, there are 1476 known interactions between 210 drugs and 204 ion channels; and in the E dataset, there are 2926 known interactions between 445 drugs and 664 enzymes. The known interactions between n drugs and m proteins are recorded by an $n \times m$ DTI matrix Z. If the ith drug is approved to target the jth protein, $Z_{i,j}=1$; otherwise $Z_{i,j}=0$. The information of the four datasets are shown in Table 1.

Table 1 The information of the benchmark datasets

The structural similarities between drugs were calculated using SIMCOMP [48] according to the size of the common substructures between two drugs.

An $n\times n$ matrix $S^d$ is used to record the similarity information between n drugs. The sequence similarity of the target proteins used the normalized version of the Smith-Waterman score [49]. Let $p_1$ and $p_2$ represent two proteins and SW(., .) be the original Smith-Waterman alignment score. The Smith-Waterman score of the normalized version of $p_1$ and $p_2$ is $s({p_1},{p_2}) = \frac{{SW({p_1},{p_2})}}{{\sqrt{SW({p_1},{p_1})} \sqrt{SW({p_2},{p_2})} }}$. An $m\times m$ matrix $S^t$ is used to store the similarity information between m target proteins.

Methods

After construction of the DTI matrix Z, the drug similarity matrix $S^d$ and a target similarity matrix $S^t$, DTI prediction is transformed into non-negative factorization of the DTI matrix with graph dual regularization terms and a prior knowledge consistency constraint. The graph dual regularization terms are used to integrate the information from the drug similarity matrix and the target similarity matrix, in order to take the intrinsic geometrical structures of the related manifolds into account. Finally, an alternating direction algorithm is used to solve the matrix factorization. Furthermore, we prove the convergence of the algorithm.

Non-negative matrix factorization

In order to obtain low-dimensional feature representations of drugs and targets in the drug–target interaction space, factorization of the DTI matrix is widely adopted. The general form of the matrix factorization is as follows:

$$\begin{aligned} Z\approx X Y^T, \end{aligned}$$

(1)

where X and Y are the latent feature matrices of drugs and targets respectively, $X \in {R^{n \times k}}, Y \in {R^{m \times k}}$, and k is the rank of feature vectors of drugs and targets ($k \ll \min (m,n)$).

To improve interpretability, the non-negativity constraint of X and Y is usually added. The optimization model of non-negative matrix factorization (NMF) is as follows:

$$\begin{array}{*{20}l} {\min \left\| {Z - XY^{T} } \right\|_{F}^{2} } \\ {s.t.~X \ge 0,Y \ge 0.} \\ \end{array}$$

(2)

Graph regularized non-negative matrix factorization

NMF aims to well approximate the DTI matrix by finding two low rank matrices, but fail to consider the geometric information in the original data. To integrate the geometric information, Cai et al. [50] proposed graph regularized non-negative matrix factorization (GNMF) which introduces a graph regularization item. The cost function of GNMF is as follows:

$$\begin{aligned} \mathcal {O}_{gr}=\frac{1}{2}\left\| {Z - X{Y^T}} \right\| _F^2&+{\lambda }{\text {Tr}}({Y^T}(Y^T(D-W)) Y)\\&s.t. ~ X \ge 0, Y \ge 0, \end{aligned}$$

(3)

where ${\text {Tr}}$ is the trace of a matrix, $\lambda$ is regularization parameter, W is the weight matrix representing a neighbor graph of the data points, and D is a diagonal matrix such that ${D_{ii}} = \sum \nolimits _l {{W_{il}}}$. The matrix $D-W$ is called graph Laplacian and denoted by $\mathcal {L}$ in the following. Furthermore, considering the geometric structure of data manifold and feature manifold, Shang et al. [43] proposed graph dual regularization non-negative matrix factorization (GDNMF), whose cost function is:

$$\begin{aligned} \begin{aligned} \mathcal {O}_{gd}&=\frac{1}{2}\left\| {Z - X{Y^T}} \right\| _F^2 +{\lambda _y}{\text {Tr}}({Y^T}(Y^T \mathcal {L}_y Y)\\&+{\lambda _x}{\text {Tr}}({X^T}(X^T \mathcal {L}_x X)\\&s.t. ~ X \ge 0, Y \ge 0. \end{aligned} \end{aligned}$$

(4)

From the similarity matrices $S^d$ and $S^t$, we could obtain the geometric information of drugs and targets. First we construct two p-nearest neighbor graphs $N^d$ and $N^t$ of drugs and targets, respectively.

For two drugs $d_i$ and $d_j$, the weight of the edge between vertices i and j in the p-nearest neighbor graph $N^d$ is defined as follows.

$$\begin{aligned} {N_{ij}^d} = \left\{ \begin{array}{l} 1, \quad j \in {\mathcal {N}_p}(i) \text { and } i \in {\mathcal {N}_p}(j) \\ 0, \quad j \notin {\mathcal {N}_p}(i) \text { and } i \notin {\mathcal {N}_p}(j) \\ 0.5, \quad \text {otherwise,} \\ \end{array} \right. \end{aligned}$$

(5)

where $\mathcal {N}_p(i)$ and $\mathcal {N}_p(j)$ are the sets of p most similar drugs of drugs $d_i$ and $d_j$ according to $S^d$, respectively. $N^d$ is used to make the drug similarity matrix $S^d$ sparse as follows.

$$\begin{aligned} {\hat{S}}_{ij}^d = {N_{ij}^d}S_{ij}^d, \forall i, j. \end{aligned}$$

(6)

${\hat{S}}^d$ is used as the weight matrix representing the drug neighbor graph. The graph Laplacian of ${\hat{S}}^d$ is ${\mathcal {L}_d} = {D^d} - {{\hat{S}}^d}$, where $D^d$ is a diagonal degree matrix with $D_{ii}^d = \sum \limits _r {{\hat{S}}_{ir}^d}$.

The same processing is performed on the target similarity matrix $S^t$ and we calculated out ${\hat{S}}^t$ the weight matrix representing the target neighbor graph as follows.

$$\begin{aligned} {\hat{S}}_{ij}^t = {N_{ij}^t}S_{ij}^t, \forall i, j. \end{aligned}$$

(7)

The graph Laplacian of ${\hat{S}}^t$ is ${\mathcal {L}_t} = {D^t} - {{\hat{S}}^t}$, where $D^t$ is diagonal degree matrix, $D_{jj}^t = \sum \limits _q {{\hat{S}}_{jq}^t}$.

Since the normalized graph Laplacian usually performs better in many actual applications, we adopted the following normalized graph Laplacian forms of $\mathcal {L}_d$ and $\mathcal {L}_t$.

$$\begin{aligned}{} & {} \begin{aligned} {\widetilde{{\mathcal {L}}}_d} = {\left( {{D^d}} \right) ^{ - 1/2}}{{{\mathcal {L}}}_d}{\left( {{D^d}} \right) ^{ - 1/2}}, \end{aligned} \end{aligned}$$

(8)

$$\begin{aligned}{} & {} \begin{aligned} {\widetilde{{\mathcal {L}}}_t} = {\left( {{D^t}} \right) ^{ - 1/2}}{{{\mathcal {L}}}_t}{\left( {{D^t}} \right) ^{ - 1/2}}. \end{aligned} \end{aligned}$$

(9)

The graph dual regularization non-negative matrix factorization (GDNMF) for DTI prediction problem is formulated as follows.

$$\begin{aligned} \begin{aligned} \mathop {\min }\limits _{(X,Y,Z)}&\frac{1}{2}\left\| {Z - X{Y^T}} \right\| _F^2 +{\lambda _d}{\text {Tr}}({X^T}\widetilde{\mathcal {L}}_d X)\\&+ {\lambda _t}{\text {Tr}}({Y^T}\widetilde{\mathcal {L}}_t Y).\\ s.t. ~&X \ge 0, Y \ge 0, \\ \end{aligned} \end{aligned}$$

(10)

where $\lambda _d$ and $\lambda _t$ are regularization parameters.

Graph dual regularized non-negative matrix factorization with prior knowledge constraint for DTI prediction

To ensure the matrix decomposition result is consistent with the prior knowledge of known DTIs, we introduce a Prior knowledge Constraint in GDNMF and formulate the DTI prediction problem as the following optimization problem (abbreviated as GRMFC).

$$\begin{aligned} \begin{aligned} \mathop {\min }\limits _{(X,Y,Z)}&\frac{1}{2}\left\| {Z - X{Y^T}} \right\| _F^2 +{\lambda _d}{\text {Tr}}({X^T}\widetilde{\mathcal {L}}_d X)\\&+ {\lambda _t}{\text {Tr}}({Y^T}\widetilde{\mathcal {L}}_t Y).\\ s.t. ~&X \ge 0, Y \ge 0, \\&{\mathrm{\mathscr {P}}_\Omega }(Z - XY^T) = 0,\\ \end{aligned} \end{aligned}$$

(11)

where $\Omega$ indexes the known drug–target interactions, i.e. the elements whose values are 1 in Z. ${\mathrm{\mathscr {P}}_\Omega }(S)$ returns a copy of S that zeros out the elements whose indices are not in $\Omega$, which is defined as follows.

$$\begin{aligned} {{\mathscr {P}_\Omega }(S)_{ij}} = \left\{ \begin{array}{l} {S_{ij}},(i,j) \in \Omega , \\ 0,(i,j) \notin \Omega . \\ \end{array} \right. \end{aligned}$$

Since the elements whose values are 1 in Z represent the known validated drug–target interactions, and we introduce the prior knowledge consistency constraint ${\mathrm{\mathscr {P}}_\Omega }(Z - XY^T) = 0$ to ensure that the matrix fraction $XY^T$ remains the known DTIs and does not lose the prior knowledge.

GRMFC is a non-convex optimization problem, and it is difficult to obtain its accurate solution. Inspired by the iteration algorithm in [41], we used an adapted alternating direction algorithm to obtain a local optimal solution. The alternating direction algorithm (ADA) is an iteration algorithm that alternatively updates X and Y. In order to use ADA to solve GRMFC efficiently, we introduce auxiliary matrices M, U and V, and transform (11) into the following equivalent form.

$$\begin{aligned} \begin{aligned} \mathop {\min }\limits _{(U,V,X,Y,M)}&\frac{1}{2}\left\| {M - X{Y^T}} \right\| _F^2 +{\lambda _d}{\text {Tr}}({X^T}\widetilde{\mathcal {L}}_d X)\\&+ {\lambda _t}{\text {Tr}}({Y^T}\widetilde{\mathcal {L}}_t Y).\\ s.t. ~&X = U, {Y^T} = V, \\&U \ge 0, V \ge 0, \\&{\mathrm{\mathscr {P}}_\Omega }(Z - M) = 0,\\ \end{aligned} \end{aligned}$$

(12)

where $U \in {R^{n \times k}}$, $V \in {R^{k \times m}}$. The auxiliary matrix M is regarded as the predicted drug–target interaction matrix.

The augmented Lagrangian of (12) is:

$$\begin{aligned} &{\mathscr {L}}(X,Y,M,U,V,\Lambda ,\Pi )\\&= \frac{1}{2}\left\| {M - X{Y^T}} \right\| _F^2+ {\lambda _d}{\text {Tr}}({X^T}{\widetilde{\mathcal {L}}}_d X) \\&+{\lambda _t}{\text {Tr}}({Y^T}{\widetilde{\mathcal {L}}}_t Y)+\Lambda \bullet (X - U)\\&+ \Pi \bullet (Y^T - V)+ \frac{\alpha }{2}\left\| {X - U} \right\| _F^2\\&+\frac{\beta }{2}\left\| {Y^T - V} \right\| _F^2, \end{aligned}$$

(13)

where $\Lambda$ and $\Pi$ are Lagrangian multipliers, $\Lambda \in {R^{n \times k}}$, $\Pi \in {R^{k \times m}}$, and $\alpha , \beta > 0$ are penalty parameters.

The alternating direction algorithm successively updates the values of matrices X, Y, M, U and V one at a time, such that $\mathscr {L}$ reaches the minimum with respect to the updated matrix while the other matrices take their most recent values. The updating rules are as follows.

$$\begin{aligned} X_{{i + 1}} & = \mathop {\arg \min }\limits_{X} {\mathscr{L}}(X,Y_{i} ,M_{i} ,U_{i} ,V_{i} ,\Lambda _{i} ,\Pi _{i} ), \\ Y_{{i + 1}} & = \mathop {\arg \min }\limits_{Y} {\mathscr{L}}(X_{{i + 1}} ,Y,M_{i} ,U_{i} ,V_{i} ,\Lambda _{i} ,\Pi _{i} ), \\ M_{{i + 1}} & = \mathop {\arg \min }\limits_{{{\mathscr{P}}_{\Omega } (Z - M) = 0}} {\mathscr{L}}(X_{{i + 1}} ,Y_{{i + 1}} ,M,U_{i} ,V_{i} ,\Lambda _{i} ,\Pi _{i} ), \\ U_{{i + 1}} & = \mathop {\arg \min }\limits_{{U \ge 0}} {\mathscr{L}}(X_{{i + 1}} ,Y_{{i + 1}} ,M_{{i + 1}} ,U,V_{i} ,\Lambda _{i} ,\Pi _{i} ), \\ V_{{i + 1}} & = \mathop {\arg \min }\limits_{{V \ge 0}} {\mathscr{L}}(X_{{i + 1}} ,Y_{{i + 1}} ,M_{{i + 1}} ,U_{{i + 1}} ,V,\Lambda _{i} ,\Pi _{i} ), \\ \Lambda _{{i + 1}} & = \Lambda _{i} + \gamma \alpha (X_{{i + 1}} - U_{{i + 1}} ), \\ \Pi _{{i + 1}} & = \Pi _{i} + \gamma \beta (Y_{{i + 1}} - V_{{i + 1}} ). \\ \end{aligned}$$

In closed form, the updating rules are as follows.

$$\begin{aligned} &{X_{i + 1}} = ({M_i}{Y_i} - {\lambda _d}{\widetilde{{\mathcal {L}}}_d}X + \alpha {U_i}- {\Lambda _i}) \\&\quad \quad {(Y_i^T{Y_i} + \alpha I)^{ - 1}},\\&{Y_{i + 1}} = ({M_{i}^T}{X_{i+1}} - {\lambda _t}{\widetilde{{\mathcal {L}}}_t}Y + \beta {V_i} - {\Pi _i})\\&\quad \quad {(X_{i+1}^T{X_{i+1}} + \beta I)^{ - 1}},\\&{M_{i + 1}} = X_{i + 1}Y_{i + 1}^T+{\mathscr {P}}_{\Omega }(M-X_{i + 1}Y_{i + 1}^T),\\&{U_{i + 1}} ={\mathscr {P}}_{+}(X_{i + 1}+{\Lambda _i}/ \alpha ),\\&{V_{i + 1}} = {\mathscr {P}}_{+}(Y_{i + 1}+{\Pi _i}/ \beta ),\\&{\Lambda _{i + 1}} = {\Lambda _i} + \gamma \alpha ({X_{i + 1}} - {U_{i + 1}}), \\&{\Pi _{i + 1}} = {\Pi _i} + \gamma \beta ({Y_{i + 1}} - {V_{i + 1}}), \end{aligned}$$

(14)

where $(\mathscr {P}_{+}(S))_{ij}=max\{S_{ij}, 0\}$, and $\gamma$ is a step length parameter which is set as 1.618 in the following experiments according to [41].

The iteration process will terminated when the changes of M smaller than a given smaller threshold. The pseudocode of the algorithm (ADA-GRMFC) is shown in Algorithm 1. The proof of convergence of ADA-GRMFC is shown in Additional file 1: Appendix.

Experiments

The performance of ADA-GRMFC has been tested in the following three aspects. First, ADA-GRMFC was compared with other state-of-the-art algorithms using the four datasets: NR, GPCR, IC and E. Second, we analyzed the effect of graph dual regularization terms on the prediction accuracy of ADA-GRMFC using ablation experiments. Third, we tested ADA-GRMFC using case studies.

To compare performances of DTI prediction algorithms, 5 repetitions of 10-fold cross-validation were performed. The final test results were the averages of the 5 repetitions of 10-fold cross-validation. The cross-validation experiments are conducted under the following two scenarios [9].

• 1. $CV_d$: The drugs are divided in ten subsets (folds), each fold is selected in turn as the test dataset and the other remained 9 folds are used as the training dataset. If the ith drug is in the test dataset, the elements in the ith row of Z are all set 0, which means the known interactions with tested drugs are removed from the input DTI matrix. It aims to evaluate the targeted protein prediction performance for the drugs without any known interactive targets.

• 2. $CV_t$: The targets are divided in ten subsets (folds), each fold is selected in turn as the test dataset and the other remained 9 folds are used as the training dataset. If the jth target in the test dataset, the elements in the jth column of Z are all set 0, which means the known interactions with tested targets are removed from the input DTI matrix. It aims to evaluate the targeting drug prediction performance for the targets without any known interactive drugs.

The area under receiver operating characteristic curve (AUC) and area under the precision-recall curve (AUPR) are used as performance evaluation metrics in the following experiments.

Comparison with state-of-the-art methods

To demonstrate the effectiveness of ADA-GRMFC in predicting DTIs, we compared ADA-GRMFC with the following seven methods, namely BLM-NII [11], WKNKN [10], RLS-WNN [8], GRMF [29], WGRMF, CMF [28], SRCMF [34], where WGRMF is a weighted form of GRMF.

Parameter settings

For BLM-NII, the combination weight $\alpha =0.5$. For WKNKN, the parameters $K=5$, $\eta =0.7$. According to the original literature of GRMF, WGRMF, CMF and SRCMF, some parameters are automated chosen using grid search [51] based on the AUPR value. Based on a previous research [29], rank k of the matrices X and Y was selected from $\left\{ {50, 100} \right\}$. For GRMF, WGRMF, CMF and SRCMF, the regularization parameter $\lambda _l$ was selected from $\{2^{-2}, 2^{-1}, 2^{0}, 2^{1}\}$. For GRMF, WGRMF and ADA-GRMFC, $\lambda _d$ and $\lambda _t$ were selected from $\{0, 10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}\}$. For CMF and SRCMF, $\lambda _d$ and $\lambda _t$ were selected from {$2^{-3}, 2^{-2}, 2^{-1}, 2^{0}, 2^{1}, 2^{2}, 2^{3}, 2^{4}, 2^{5}$}.

In terms of ADA-GRMFC, the optimal parameters combination of $\alpha$, $\beta$ and $\gamma$ are different under different experiment scenario, which makes the parameter settings more complicated. When tuning the parameter $\alpha$, we set other parameters as their optimal values. The same settings were applied to $\beta$ and $\gamma$. Figures 1, 2 and 3 in the appendix file showed the impact of the parameter $\alpha$, $\beta$ and $\gamma$ on AUC and AUPR, respectively. When the parameters change, the performances of ADA-GRMFC varied more significantly on NR datasets and GPCR datasets than on IC datasets and E datasets. The results show that the best performance was achieved when $\alpha = 0.5$, $\beta =0.01$ and $\gamma =1.618$. Thus, in the following experiments, for ADA-GRMFC, $\alpha =0.5$, $\beta =0.01$, $\gamma = 1.618$. The parameters $\epsilon =10^{-6}$.

Prediction results

Under the $CV_d$ scenario, ADA-GRMFC performs better than other methods in terms of AUC and AUPR on NR, IC, and E datasets. The AUC values of ADA-GRMFC are 0.860748, 0.798762, and 0.834382 on NR, IC, and E datasets, respectively. The AUPR values of ADA-GRMFC are 0.574956, 0.374033, and 0.39878 on NR, IC, and E datasets, respectively. On the GPCR dataset, WGRMF achieve the highest AUC and AUPR values, which are 0.868548 and 0.410652, respectively. The AUC and AUPR values of ADA-GRMFC are higher than those of other algorithms except WGRMF. The AUC and AUPR values of the different algorithms on the four datasets are shown in Tables 2 and 3, respectively. The AUC and AUPR histograms with error bars of different algorithms are shown in Figure 1a, b, respectively. The receiver operating characteristic (ROC) curves and the precision-recall (PR) curves of different methods on the four datasets are shown in Figs. 2 and 3, respectively.

Table 2 AUC values of different algorithms under $CV_d$ scenario

Table 3 AUPR values of different algorithms under $CV_d$ scenario

Fig. 1

AUC values and AUPR values of the methods on the four datasets under $CV_d$. a Histogram with error bars of AUC. b Histogram with error bars of AUPR

Fig. 2

ROC curves for different methods are plotted together under CVd on NR dataset, GPCR dataset, IC dataset, E dataset, respectively

Fig. 3

PR curves for different methods are plotted together under $CV_d$ on NR dataset, GPCR dataset, IC dataset, E dataset, respectively

Under the $CV_t$ scenario, the AUC and AUPR values of ADA-GRMFC are higher than the other methods on the four datasets. The AUC values of ADA-GRMFC are 0.799721, 0.896419, 0.948086, and 0.939765 on NR, GPCR, IC, and E datasets, respectively. The AUPR values of ADA-GRMFC on NR, GPCR, IC, and E datasets are 0.454528, 0.598742, 0.812833, and 0.806995, respectively. The AUC values and AUPR values of different algorithms on the four datasets are shown in Tables 4 and 5, respectively. The AUC and AUPR histograms with error bars of different algorithms are shown in Fig. 4a, b, respectively. ROC and PR curves of different algorithms are shown in Figs. 5 and 6 on the four datasets, respectively.

Table 4 AUC values of different algorithms under $CV_t$ scenario

Table 5 AUPR values of different algorithms under $CV_t$ scenario

Fig. 4

AUC values and AUPR values of the methods on the four datasets under $CV_t$. a Histogram with error bars of AUC. b Histogram with error bars of AUPR

Fig. 5

ROC curves for different methods are plotted together under $CV_t$ on NR dataset, GPCR dataset, IC dataset, E dataset, respectively

Fig. 6

PR curves for different methods are plotted together under $CV_t$ on NR dataset, GPCR dataset, IC dataset, E dataset, respectively

Ablation experiments

Our model includes two graph dual regularization terms: a regularization term for drugs and regularization term for targets. To evaluate the impact of the graph dual regularization terms on the performance of ADA-GRMFC, we conduct ablation experiments on the benchmark datasets.

In Tables 6, 7, 8 and 9, parameter $\lambda _d$ and $\lambda _t$ of ADA-GRMFC are chosen using grid search. ADA-GRMFC with $\lambda _d=0$ means that the graph regularization term for drugs is not used. ADA-GRMFC with $\lambda _t=0$ means that the graph regularization terms for targets is not used. When we use regularization terms for drugs and targets, ADA-GRMFC has the highest prediction performance in $CV_d$ and $CV_t$. In $CV_d$, when $\lambda _d=0$, the values of AUC and AUPR of ADA-GRMFC are significantly decreased. The AUC values have decreased by 23%, 19%, 44%, 36% on NR, GPCR, IC and E datasets, respectively. The AUPR values have decreased by 92%, 87%, 91%, 97% on NR, GPCR, IC and E datasets, respectively. Similarly, in $CV_t$, if the graph regularization terms for targets is not used, the performances of ADA-GRMFC is significantly decreased too. When $\lambda _t=0$, the AUC values have decreased by 34%, 33%, 43%, 45% on NR, GPCR, IC and E datasets, respectively. The AUPR values have decreased by 77%, 91%, 95%, 98% on NR, GPCR, IC and E datasets, respectively. These results show that regularization terms for drugs and targets contribute the improvement of DTI prediction performance of ADA-GRMFC significantly.

Table 6 AUC results for ADA-GRMFC variants under $CV_d$

Table 7 AUPR results for ADA-GRMFC variants under $CV_d$

Table 8 AUC results for ADA-GRMFC variants under $CV_t$

Table 9 AUPR results for ADA-GRMFC variants under $CV_t$

Case studies

To further evaluate the ability of ADA-GRMFC to find new targets for a drug and new drugs for a target in practice, two case studies concerning the drug olanzapine and the target estrogen receptor alpha were conducted.

In the first case study, we predicted targets that interact with the drug olanzapine on the G protein-coupled receptors (GPCR) dataset using ADA-GRMFC. Olanzapine is an antipsychotic drug which could target many receptors, and it was recently found that olanzapine could be an attractive antiemetic drug [52]. The known interactions of olanzapine with targets were deleted from the the training dataset, and the candidate targets of olanzapine predicted by ADA-GRMFC were prioritized according to the prediction scores. At last, the top 10 highest-scoring predicted targets were picked out to be validated using the databases KEGG [53] and DrugBank [54]. The results showed that all 10 targets were correctly predicted. The detailed results of the predictions are shown in Table 10.

In the second case study, we predicted candidate drugs for the target estrogen receptor alpha (ER$\alpha$) on the NR dataset and aimed to assess the ability of ADA-GRMFC to predict candidate drugs for targets with no known targeting drugs. ER$\alpha$ is mainly expressed in reproductive tissues (uterus, ovaries), breast, kidney, bone, white adipose tissue and liver, and is over-expressed in more than half of breast cancers [55]. The known interactions of ER$\alpha$ with drugs were removed from the training dataset, and the candidate drugs of ER$\alpha$ predicted by ADA-GRMFC were prioritized according to the prediction scores. The top 20 highest-scoring predicted drugs were selected to be validated against the databases KEGG and DrugBank. Among the predicted 20 drugs, 17 drugs had evidences to target ER$\alpha$. The detailed results of the case study are shown in Table 11.

Table 10 Top 10 predicted targets of olanzapine by ADA-GRMFC on the NR dataset

Table 11 Top 20 predicted drugs targeting estrogen receptor alpha by ADA-GRMFC on the GPCR dataset

Conclusion

The knowledge of interactions between drugs and targets could help to find the novel usage of drugs. In the paper, we propose a matrix factorization based method, ADA-GRMFC, to predict interactions between drugs and targets. ADA-GRMFC uses graph dual regularization terms to capture structural information from the drug similarity matrix and the target similarity matrix. At the same time, the prior knowledge consistency constraint is used to ensure the matrix decomposition result is consistent with the known DTIs. Finally, an alternating direction algorithm is used to solve the matrix factorization with constraints. Extensive experiments show that ADA-GRMFC outperforms the state-of-the-art methods in predicting DTIs.

The reasons for the superior performance of ADA-GRMFC are as follows. First, unlike traditional matrix factorization algorithms, the prior knowledge consistency constraint ensures that the matrix decomposition result is consistent with the prior knowledge of known DTIs. Second, the graph dual Laplace regular terms not only overcome overfitting of model, but also obtain underlying structural information about the data. Finally, we use alternating direction algorithm with fast convergence to solve the constrained problem.

However, ADA-GRMFC also has limitations. The lack of known drug–target associations may affect the performance of ADA-GRMFC, and including more information related with drugs and targets would help to improve the prediction ability. The values of the parameters of ADA-GRMFC are set by grid search which is time consuming, and appropriate methods to choose optimal parameters need further research.