In this work, a new computational model called GNMFDMA is developed to predict associations of drug with miRNA. The GNMFDMA approach can be summarized into the following three steps (See Figure 1). First, the similarity matrix of drugs is constructed according to the drug chemical structure similarity, indication phenotype similarity of drug, drug side effect similarity and gene functional consistency similarity of drug. The similarity matrix of miRNAs is constructed based on gene functional consistency and disease indication phenotype similarity of miRNA. Second, to extend GNMFDMA to novel drugs and miRNAs, we use weighted K nearest neighbor profiles to re-construct the drug-miRNA association adjacency matrix. Finally, graph Laplacian regularization collaborative standard non-negative matrix factorization is utilized to discover drug-miRNA potential associations.

In order to infer potential associations of drug with miRNA using non-negative matrix factorization, we construct the drug-drug interaction network and miRNA-miRNA interaction network by integrating the four categories of drug-drug similarities and two categories of miRNA-miRNA similarities, respectively. Besides, drug-miRNA association network is constructed using the known drug-miRNA association pairs.

In this work, the verified 664 drug-miRNA associations were obtained from the SM2miR database, which can be accessible at http://bioinfo.hrbmu.edu.cn/SM2miR/ (Liu et al., 2012). In the 664 known associations, 831 drugs obtain from PubChem (Wang et al., 2009), DrugBank (Knox et al., 2010), and SM2miR; 541 miRNAs collect from PhenomiR (Ruepp et al., 2010), HMDD (Lu et al., 2008), miR2Disisease (Jiang et al., 2008), and SM2miR databases. However, there are some drugs and miRNAs without any known association information. For this reason, these drugs and miRNAs are deleted, and the duplicated entries are also removed. After screening, the drug-miRNA association network with 664 different associations is constructed for prediction, including 39 drugs and 286 miRNAs (See Table 1). Based on the drug-miRNA association network, the original association adjacency matrix Y ∈ R n × m is defined, where m and n represent the number of miRNAs and drugs, respectively. The element value Y i , j is set to one if drug d i is confirmed to be associated with miRNA m j , otherwise it is 0.

Previous studies have shown that similarities based on chemical structure (Hattori et al., 2003), indication phenotype (Gottlieb et al., 2011), side effect (Gottlieb et al., 2011) and gene functional consistency (Lv et al., 2011) are effectively tools to infer the relationships between drugs. In this work, to avoid the bias of single similarity measurement and contribute the discovery of new interactions, four types of drug similarity were integrated according to the model of Lv et al. (2015). The four types of drug similarity are drug chemical structure similarity, disease indication phenotype similarity of drug, drug side effect similarity and gene functional consistency similarity of drug, respectively (Lv et al., 2015). We use matrix S g d to represent the drug similarity information based on gene functional consistency. The element S g d i , j of matrix S g d is the functional consistency similarity of drug d i and drug d j . At the same time, S c d , S s d and S d d denote the similarity matrices based on chemical structure, side effect and disease indication phenotype, respectively. For each pair of drugs, four types of similarity are combined to calculate the overall similarity as follows: S D = ω 1 S c d + ω 2 S d d + ω 3 S g d + ω 4 S s d ω 1 + ω 2 + ω 3 + ω 4 (1) where the weight value ω 1 , ω 2 , ω 3 and ω 4 are assigned as 1, respectively. The size of S D is n × n , The element S D i , j denotes the similarity of drug d i with drug d j .

The similarity of miRNA is constructed in this work using the model proposed by Lv et al. (2015), which is based on disease indication phenotype similarity of miRNA and gene functional consistency similarity of miRNA, respectively (Gottlieb et al., 2011; Lv et al., 2011). S d m and S g m denote disease indication phenotype similarity of miRNA and gene functional consistency similarity of miRNA. Then, we calculate the overall similarity of miRNA by integrating the two types of similarity S g m and S d m as follows: S M = σ 1 S d m + σ 2 S g m σ 1 + σ 2 (2) where the weight value σ 1 and σ 2 are assigned as 1, respectively. The size of S M is m × m , the element S M i , j is the similarity of miRNA m i with miRNA m j .

Let D = d 1 , d 2 , ⋯ , d n and M = m 1 , m 2 , ⋯ , m m are the set of n drugs and m miRNAs. The i t h row vector Y d i = Y i 1 , Y i 2 , ⋯ , Y i m and the j t h column vector Y m j = Y 1 j , Y 2 j , ⋯ , Y n j of matrix Y denote the interaction profiles of drug d i and miRNA m j , respectively. For a novel drug without any known associated miRNAs or a novel miRNA without any known associated drugs, there are no interactions in their profiles. In fact, many of unknown drug-miRNA association pairs (or 0’s) in Y could be potential true associations, which may result in a higher false positive rate and reduce prediction performance. In order to address this problem, a preprocessing step (WKNKN) is performed to construct new interaction profiles based on their known neighbors.

For each drug d l , all other drugs are ranked in descending order on the basis of their similarity to d l . Then, the new interaction profile for drug d l is obtained based on their corresponding interaction profiles of the K known drugs nearest to d l (Ezzat et al., 2017): Y d d l = 1 ∑ 1 ≤ i ≤ K S D d i , d l ∑ i = 1 K θ i Y d i (3) where θ i = α i − 1 * S D d i , d l is the weight coefficient, a larger θ i represents that d i and d l are more similar. α ∈ 0,1 is a decay term. The same procedure for miRNA, for each miRNA m p , the new interaction profile can be defined as follows: Y m m p = 1 ∑ 1 ≤ j ≤ K S M m j , m p ∑ j = 1 K θ j Y m j (4) Similarly, all other miRNAs are ranked in descending order according to their similarity to m p . θ j = α j − 1 * S M m j , m p is the weight coefficient.

Then, we merge the two matrices of Y d and Y m , and replace Y i j = 0 with the associated likelihood score. Finally, the novel drug-miRNA association adjacency matrix is obtained: Y = max ⁡   Y ,   Y d m (5) where Y d m = Y d + Y m 2 (6)

Based on the spectral graph and manifold learning theories that the nearest neighbor graph can maintain the local geometry of the original data points, and the sparseness technique of similarity matrix has been successfully applied in graph regularization (Cai et al., 2010; You et al., 2010; Li et al., 2016b). At the same time, the drugs and miRNAs located in the same cluster often have more similar functions. Thus, we calculate the affinity graphs ( S D * ; S M * ) for drug space and miRNA space using p -nearest neighbor. Then, the weight matrix of drug is defined according to the drug similarity matrix S D as follows: G i j D =   1   i ∈ N p d j & j ∈ N p d i 0   i ∉ N p d j & j ∉ N p d i 0.5   o t h e r w i s e   (7) where N p d i and N p d j are the sets of p-nearest neighbors of drug d i and drug d j , respectively. Finally, the sparse similarity matrix S D * for drugs is calculated as: ∀ i , j ,   S i j D * = S i j D G i j D (8)

Similarly, the sparse similarity matrix S M * for miRNAs is calculated as follows: ∀ i , j ,   S i j M * = S i j M G i j M (9)

Non-negative matrix factorization (NMF) method has been effectively applied for data representation. NMF decomposes an original matrix into two non-negative matrices whose product is as equal to the original matrix as possible. At the same time, it can also achieve the purpose of dimensionality reduction. In this work, NMF is used to decompose the drug-miRNA association adjacency matrix Y n × m into W k × n ; H k × m ( k < min ⁡ ⁡ m ,   n ), and Y ≅ W T H . The problem of drug-miRNA association prediction can be expressed by the following objective function: min W , H Y − W T H F 2   s , t .   W ≥ 0 ,   H ≥ 0 (10) where ∙ F 2 is the Frobenius norm and k is the subspace dimensionality. However, in the Euclidean space, the intrinsic geometrical of the drug or miRNA space cannot be discovered by standard NMF (Wang et al., 2022b). To prevent overfitting and enhance generalization capability of the model, the graph Laplacian regularization terms and Frobenius norm regularization terms (Tikhonov L 2 ) are introduced to the standard NMF. The graph Laplacian regularization can ensure local invariance for data space (Cai et al., 2010). Here, we use graph Laplacian regularization to ensure close miRNAs or drugs to be adequately close to each other in miRNA or drug corresponding space. In addition, the Frobenius norm regularization terms are utilized to guarantee the smoothness of W and H . Therefore, the objective function of GNMFDMA can be transformed into: min W , H Y − W T H F 2 + λ ∑ i ≤ j = 1 n w i − w j 2 S i j D * + ∑ i ≤ j = 1 m h i − h j 2 S i j M *   + β W F 2 + H F 2         s , t .   W ≥ 0 ,   H ≥ 0 (11) where β and λ represent the sparseness constraint coefficient and regularization coefficient, respectively. w i and h j are the i t h and j t h columns of W and H , respectively. R d = ∑ i ≤ j = 1 n w i − w j 2 S i j D * = ∑ j = 1 n w j T w j ∑ i , j = 1 n S i j D * − ∑ i , j = 1 n w i T w j S i j D * = ∑ j = 1 r w j T w j D j j − ∑ i , j = 1 r w i T w j S i j D * = T r W D d W T − T r W S D * W T = T r W L d W T (12) and R m = ∑ i ≤ j = 1 m h i − h j 2 S i j M * = T r H L m H T (13) here, T r ∙ denotes the trace of matrix. R d and R m are the graph Laplacian regularization terms. L d = D d − S D * is graph Laplacian matrix of S D * , L m = D m − S M * is graph Laplacian matrix of S M * , respectively (Liu et al., 2014). D d and D m are the diagonal matrices, D d i , i = ∑ l = 1 n S i l D * and D m j , j = ∑ p = 1 m S j p M * . Eq. 11 can be expressed as follows: min W , H Y − W T H F 2 + β W F 2 + H F 2 + λ T r W L d W T + λ T r ( H L m H T ) = T r ( Y Y T ) + T r ( W T H H T W ) − 2 T r ( Y H T W ) + β T r W T W + β T r H T H + λ T r W L d W T + λ T r H L m H T (14)

To minimize Eq. 14, we introduce Lagrange multipliers method to solve this problem. Let Lagrange multipliers ψ = φ k i and Φ = ϕ k j to ensure w k i ≥ 0 and h k j ≥ 0 . The corresponding optimization function F of Eq. 14 is formularized as: F = T r ( Y Y T ) + T r ( W T H H T W ) − 2 T r ( Y H T W ) + β T r ( W T W ) + β T r H T H + λ T r W L d W T + λ T r H L m H T + ψ T r W T + Φ T r H T (15)

The partial derivatives of F for W and H are: ∂ F ∂ W = 2 H H T W − 2 H Y T + 2 β W + 2 λ W L s + ψ (16) ∂ F ∂ H = 2 W W T H − 2 W Y + 2 β H + 2 λ H L m + Φ (17)

Then, the Karush–Kuhn–Tucker (KKT) condition φ k i w k i = 0 ; ϕ k j h k j = 0 are used in Eq. 16 and Eq. 17 (Facchinei et al., 2014). We can obtain the following Equations: H H T W k i w k i − H Y T k i w k i + β W k i w k i + λ W D d − S D * k i w k i = 0 (18) W W T H k j h k j − W Y k j h k j + β H k j h k j + λ H D m − S M * k j h k j = 0 (19)

Thus, the updating rules for w k i and h k j can be obtained as follows: w k i ← w k i H Y T + λ W S D * k i H H T W + β W + λ W D d k i (20) h k j ← h k j W Y + λ H S M * k j W W T H + β H + λ H D m k j (21)

Updating w k i and h k j with Eq. 20 and Eq. 21 until W and H reach the following convergence conditions: ∀ i ,   w i l + 1 − w i l F 2 ≤ 10 − 4 (22) ∀ j ,   h j l + 1 − h j l F 2 ≤ 10 − 4 (23)

Ultimately, the predicted drug-miRNA associations adjacency matrix Y * is calculated by Y * = W T H . The elements of matrix Y * are regarded as the drug-miRNA association predicted scores. For each drug-miRNA pair, all the miRNAs are sorted in descending order based on the predicted scores. In theory, the top ranked miRNAs in predicted matrix Y * are more possible to be related to the corresponding drug.

To systematical evaluate the performance of GNMFDMA, we carry out five-fold cross-validation (5-CV) experiments on SM2miR database and compare it with four state-of-the-art predictors: SMiR-NBI (Li et al., 2016a), SLHGISMMA (Yin et al., 2019), TLHNSMMA (Qu et al., 2018) and RFSMMA (Wang et al., 2019b). Specifically, in the framework of five-fold cross-validation, 664 known drug-miRNA association pairs are randomly divided into five equal subsets. Four subsets of them are taken in turn as the training samples to train the prediction model, and the remaining one subset is regarded as the test sample. In this work, the AUC values (the area under the ROC curve) are used to assess the prediction performance of various models. AUC = 0.5 represents randomly prediction, whereas AUC = 1 represents that the prediction performance of the method is perfect.

In this paper, the parameter values are chosen by 5-CV experiment on the training dataset. GNMFDMA has the following five parameters, the neighborhood size K and decay value α are chosen from 1 ,   2 ,   3 ,   4 ,   5 and 01 ,   0.2 ,   0.3 , ⋯ ⋯ ,   0.9 ,   1 when the adjacency matrix is reformulated, respectively. For non-negative matrix factorization, three parameters are subspace dimensionality k , regularization coefficient λ and sparseness constraint coefficient β , whose combinations are regarded from the following values: k ∈ 15 ,   20 ,   25 ,   30 ,   35 , λ ∈ 0.2 ,   0.6 ,   1 ,   2 and β ∈ 0.002 ,   0.02 ,   0.2 ,   0.6 . According to previous studies (Cai et al., 2010), let p = 5 when constructing the graph spaces for drug and miRNA. In order to more fairly comparison with previous methods, the parameters in other methods are all taken the optimal values recommended by authors. Finally, the parameters optimized values of our model are K = 3 , α = 0.9 , k = 35 , λ = 1 and β = 0.02 .

The performance of GNMFDMA is evaluated by comparing with the previous computational models: SMiR-NBI, SLHGISMMA, TLHNSMMA and RFSMMA. For the above methods, we all use 5-CV to evaluate their performance. Figure 2 draws the ROC curves of GNMFDMA, Table 2 displays the AUC values of all compared approaches. The AUC values of GNMFDMA, SMiR-NBI, SLHGISMMA, TLHNSMMA and RFSMMA are 0.9193, 0.7104, 0.7724, 0.8168 and 0.8389, respectively. GNMFDMA achieves the best performance, which are 20.89%, 14.69%, 10.25% and 8.04% higher than the other four computational methods, respectively.

Additionally, in order to calculate the ratio of exact identifications in the predicted results, sensitivity (Sen), accuracy (Acc), precision (Pre) and F1-Score are widely applied to measure the model performance. S e n . = T P T P + F N (24) P r e . = T P T P + F p (25) A c c . = T N + T P T N + T P + F N + F p (26) F 1 − S c o r e = 2 × P r e . × S e n . P r e . + S e n . (27)

Here, when given a cutoff value, TP and FP denote the number of true positive samples and false positive samples, whose prediction scores higher than cutoff value; TN and FN are the number of true negative samples and false negative samples, whose prediction scores lower than cutoff value. In this work, the threshold of specificity is set 85% to calculate sensitivity, accuracy, precision and F1-Score, respectively. Table 3 exhibits the sensitivity, accuracy, precision, and F1-Score by GNMFDMA under 5-CV.

In general, the predicted results obtained from top-ranked are more convincing compared with those obtained from other portions. The more true association pairs that are correctly retrieved from the top-ranked, the predictor is more effective. For this reason, we calculate the correct recovery of association pairs at different thresholds when all 664 known drug-miRNA association pairs are used as training samples. The top 10%, 15% and 20% drug-related miRNAs in prediction result, GNMFDMA correctly retrieved 429 (64.61%), 532 (80.12%) and 617 (92.92%) association pairs, respectively. The comparison between the original association adjacency matrix and the predicted association matrix is shown in Figure 3. These results show that GNMFDMA can effectively retrieve the true association pairs with a lower false negative rate. In summary, the method of GNMFDMA has powerful ability for identifying drug-associated miRNAs.

In order to investigate the effects of preprocessing step (WKNKN) for GNMFDMA, we compared the performance of GNMFDMA and GNMFDMA* under 5-CV. For GNMFDMA, we implement a preprocessing step (WKNKN) to re-construct the drug-miRNA association adjacency matrix based on their known neighbors before performing non-negative matrix factorization, which can supplement more interaction information to give assistance for predicting new drugs and miRNAs. In addition, the preprocessing step is also helpful for predicting those drugs or miRNAs with sparse known associations. For GNMFDMA*, the preprocessing step is ignored and matrix factorization is directly performed on the original adjacency matrix for inferring drug-associated miRNAs. Figure 2 and Figure 4 represent the ROC curves of GNMFDMA and GNMFDMA* under 5-CV, the AUC values achieved by GNMFDMA and GNMFDMA* are 0.9193 and 0.8507, respectively. The results demonstrate that the performance of GNMFDMA is significantly improved after performing the preprocessing step.

To further demonstrate the availability of GNMFDMA to discover potential associations of drug with miRNA, case studies are conducted for three common small molecule drugs, 5-Aza-CdR, 5-FU and Gemcitabine based on the SM2miR dataset. In each independent case study, all known 5-Aza-CdR (5-FU or Gemcitabine)-related miRNAs are removed (all miRNAs are regarded as the potential candidates of corresponding drug), the remaining known associations are utilized as the training samples. Next, for each investigated drug, these miRNAs are sorted in descending order according to the predicted scores, it means that the top-ranked miRNAs tend to be related to the corresponding drug.

We use the experimental literature to verify the predicted potential miRNAs for three corresponding drugs. The top 50 potential candidate miRNAs associated with 5-Aza-CdR, 5-FU and Gemcitabine predicted by GNMFDMA are exhibited in Table 4, Table 5 and Table 6, respectively. 30, 31 and 34 out of the top-50 miRNAs inferred by GNMFDMA are verified to be related to the corresponding drug by the experimental literature, respectively. For example, the expression of hsa-let-7d and hsa-let-7e was significantly down-regulated in gemcitabine-resistant cells (Li et al., 2009). Up-regulation of has-let-7 by natural agents can lead to the reversal of epithelial-to-mesenchymal transition in gemcitabine-resistant pancreatic cancer cells. Hsa-miR-125a promotes chemical resistance of pancreatic cancer cells to Gemcitabine by targeting A20 (Yao et al., 2016). In addition, the SM2miR database confirmed that hsa-miR-125a is also associated with drug 5-Aza-CdR. That is, one miRNA may be targeted by multiple small molecule drugs. The above results show that GNMFDMA can effectively predict new drugs or miRNAs without any known relationships, which has important reference significance for related biomedical experiments.