Here, to predict potential associations between microbes and diseases, the model of GRNMFHMDA (See Figure 1) can be decomposed into three steps: (1) data preparation, in which adjacency matrix, microbe similarity, and disease similarity were established; (2) the preprocessing step, in which unknown microbe-disease pairs were assigned with associated likelihood scores based on the calculation of weighted K nearest neighbor profiles for microbes and diseases; (3) GRNMF, in which Tikhonov (L₂) and Graph Laplacian regularization were introduced into the standard NMF framework to obtain the final score matrix.

From the Human Microbe-Disease Association Database (HMDAD, http://www.cuilab.cn/hmdad) (Ma et al., 2017), we can download 483 known microbe-disease associations between 292 microbes and 39 human diseases. However, since some microbe-disease associations we downloaded are the same, there were only 450 known associations after removing the duplicate parts according to different evidences. In order to represent the associations information in a more convenient and efficient way, we defined an adjacency matrix Y ∈ R^m*n, where m and n denoted the number of microbes and diseases, respectively. Moreover, the element Y(m_i, d_j) was set to 1 if microbe m_i and disease d_j had known association, otherwise 0.

There is a hypothesis that similar microbes (i.e., microbes exhibiting a similar pattern of interaction and non-interaction with the diseases of a microbe-disease association network) are inclined to be associated with the same disease, on which many previous studies had relied to construct the Gaussian interaction profile kernel similarity for microbes (Chen et al., 2017a; Huang Z. A. et al., 2017). In this article, based on the same assumption, we first represent the interaction profile for each microbe with a binary vector involving the association information between the microbe and each disease in the known microbe-disease association network. On the basis of the definition of adjacency matrix Y, the ith row vector (Y(m_i) = (Y_i1, Y_i2, …, Y_in)) can be used to denote the interaction profile of microbe m_i. Thus, according to the method of van Laarhoven et al. (2011), the Gaussian interaction profile kernel similarity between microbe m_i and m_j can be defined as follows:

where

Here, γ_m is the adjustment coefficient that can be obtained by normalizing another bandwidth parameter γ′_m.

The construction of the Gaussian interaction profile kernel similarity for diseases is based on the assumption that similar diseases (i.e., diseases exhibiting a similar pattern of interaction and non-interaction with the microbes of a microbe-disease association network) are more likely to be associated with similar microbes. Here, the interaction profile for each disease is also represented by a binary vector containing the association information between the disease and each microbe in the known microbe-disease association network. Based on the same method of van Laarhoven et al. (2011), the jth column vector (Y(d_j) = (Y_1j, Y_2j, …, Y_mj)) denotes the interaction profile of disease d_j and the Gaussian interaction profile kernel similarity between disease d_i and d_j can be defined as follows:

where

Similarly, γ_d is the adjustment coefficient that can be calculated by normalizing another bandwidth parameter γ′_d.

As we have mentioned above, Gaussian interaction profile kernel similarity is used in our model to measure the similarity of microbes and diseases. However, since the Gaussian interaction profile kernel similarity is an association information-based measurement, it is essential to combine more types of microbe or disease similarities based on other available biological information. Indeed, according to different biological data, many researchers have developed their own method to measure the similarity of microbes or diseases. For instance, Zhou et al. (2014) proposed a model of symptom-based human disease network (HSDN) to measure the disease similarity based on co-occurrence of disease/symptom terms recorded in different literatures. In this work, we implemented HSDN to calculate symptom-based disease similarity (SDM) and then constructed a new disease similarity matrix (S^d) by integrating SDM with S^d′ in an average way according to the study of Chen et al. (2017a):

Due to the fact that values in interaction profiles of microbes or diseases without known associations are all zeros, the prediction performance may be affected to some extent. Considering that, to deal with the above mentioned problem, we came up with a preprocessing step to establish new interaction profiles both for microbes and diseases. For each microbe m_q, we first find out its K nearest known microbes, each of which must has at least one known association. Next, the similarity information between m_q and its K nearest known microbes together with the information of their corresponding K interaction profiles are combined to calculate the new interaction profile as follows:

where

Here, m₁ to m_K denote the K nearest known microbes of m_q which were sorted in descending order based on the similarity values between them. The function of the weight coefficient w_i is that the corresponding similarity value is assigned higher weight if m_i is more similar to m_q. Besides, α is a decay term whose value is in the range of [0,1] and Q_m is the normalization term.

In a similar way, the new interaction profile for each disease d_p can be defined as follows:

After calculating the new interaction profiles from microbe perspective and disease perspective, we combine Y_m and Y_d as follows:

where a₁ and a₂ are two weight coefficient and both of them are set to 1 for simplicity.

Finally, to replace the element Y(m_i, d_j) = 0 with an associated likelihood score, we use the following equation to update the original adjacency matrix Y.

As a common method, the purpose of the standard NMF is to find two non-negative matrices whose product is an optimal approximation to the original matrix (Sotiras et al., 2015; Xu et al., 2015). Therefore, the adjacency matrix Y ∈ R^m*n can be decomposed into two parts after implementing NMF, namely, W ∈ R^m*k and H ∈ R^n*k (Y ≈ WH^T). Accordingly, we can further get the following standard optimization problem:

where L(W, H) is a regularization term to prevent overfitting.

Here, motivated by the study of Xiao et al. (2017) and the standard NMF framework, we introduced other two terms, the Tikhonov (L₂) (Guan et al., 2011) and graph Laplacian regularization (Cai et al., 2011), to predict microbe-disease associations. The utilizing of Tikhonov regularization aims to obtain a smooth solution (W and H), while the purpose of introducing graph regularization is to ensure a part-based representation through taking full advantage of the data geometric structure. Thus, we can construct the optimization problem of GRNMF as follows:

Here, λ_l, λ_m and λ_d are the corresponding regularization coefficients. Besides,w_i and h_j are defined as ith rows of W and jth rows of H, respectively. In order to avoid negative affects to the prediction performance of our model, we introduced sparse weight matrices of S^d* and S^m*that are constructed on the basis of the geometrical information of disease and microbe data spaces (S^d and S^m), respectively. Then, Equation (14) can be transformed into:

where Tr(•) represents the trace of a matrix. Here, L_m and L_d are the corresponding graph Laplacian matrices for S^m* and S^d* that can be calculated as follows:

where D_m and D_d are the diagonal matrices whose entries are row (or column) sums of S^m* and S^d*, respectively.

Based on the information of the nearest neighbor graph on a scatter of data points, researchers came up with a conclusion that local geometric structure is able to be effectively modeled (Cai et al., 2011; Li et al., 2017). Since microbes or diseases appearing in the same cluster are more likely to behave similarly, according to the above conclusion, we construct the graph matrices S^m* and S^d* in terms of microbe space and disease space respectively on the basis of the p nearest neighbors and corresponding clustering information. Here, we use the ClusterONE method (Nepusz et al., 2012) to construct the graph S^m*from microbe space, in which the weight matrix X^m is generated based on the microbe similarity matrix S^m as follows:

where N(m_i) and N(m_j) are the sets of p nearest neighbors of m_i and m_j, respectively. C denotes to any one of the clusters obtained by ClusterONE method and we define the graph matrix S^m* for microbes as follows:

In a similar way as the computation of S^m*, we calculate the graph matrix S^d* according to the disease similarity matrix S^d.

Here, we defined Φ = [ φ_ik ] and Ψ = [ ψ_jk ] as the Lagrange multipliers for the constrains w_ik ≥ 0 and h_jk ≥ 0, respectively. In this work, we first convert the optimization problem in Equation (15) to an unconstraint problem, then minimize this problem by utilizing the corresponding Lagrange function L_f as follows:

To solve the above problem, we first calculate the partial derivatives with respect to W and H as follows:

After using the Karush-Kuhn-Tucker (KKT) conditions of φ_ikw_ik = 0 and ψ_jkh_jk = 0 (Facchinei et al., 2014), we can obtain the equations for w_ik and h_jk as follows:

Finally, on the basis of the above two equations, we can get the updating rules for w_ik and h_jk as follows:

Based on the above two updating formulas, we can obtain the final two non-negative matrices W and H until convergence. Subsequently, we calculate the score matrix Y^* for microbe-disease pairs by utilizing Y^* = WH^T, in which the higher score of a microbe-disease pair indicates that the microbe is more likely to be associated with the corresponding disease. In addition, for better understanding, we provided the pseudocode of the whole GRNMF algorithm (See Figure 2).

Cross validation, a widely used assessment method, was introduced to evaluate the prediction performance of GRNMFHMDA. In this study, we utilized two types of cross validations, namely, global LOOCV and local LOOCV. For the global LOOCV, each of the known microbe-disease associations was in turn considered to be the test sample while the remaining known associations were treated as the training samples. Besides, all of the unknown microbe-disease pairs were regarded as the candidate samples which would be used in the ranking process. After implementing GRNMFHMDA, we ranked each test sample with all candidate samples according to their predicted scores. As for local LOOCV, the difference is that the test sample was only ranked with the candidate samples involving the investigated disease.

In each cross validation process, we would consider that the test sample was successfully predicted if the ranking of the test sample was higher than the given threshold. Further, based on the ranks of all test samples, we drew a receiver operating characteristic (ROC) curve through calculating the ratio between true positive rate (TPR, sensitivity) and false positive rate (FPR, 1-specificity) under different thresholds both for global LOOCV and local LOOCV. Sensitivity meant the ratio between the number of test samples ranking higher than the given threshold and the number of positive samples (known microbe-disease associations), while 1-specificity denoted the percentage of the number of negative microbe-disease pairs whose ranks were lower than the given threshold. Moreover, area under the ROC curve (AUC) was calculated to make quantitative evaluation for our model's prediction performance. The model would be considered to be able to perfectly predict all associations if the value of AUC equaled to 1, while the model was only supposed to be able to make random prediction if the value of AUC equaled to 0.5. As a result, GRNMFHMDA obtained AUCs of 0.8715 and 0.7898 in global LOOCV and local LOOCV, respectively. Furthermore, the prediction performance of our model outperformed the KATZHMDA both in global LOOCV (0.8644) and local LOOCV (0.6998), which proved the superior accuracy and reliability of our model in predicting microbe-disease associations (See Figure 3).

Here, we put forward two types of case studies on three different common human diseases with the purpose of further assessing the prediction performance of GRNMFHMDA. On the basis of the known microbe-disease associations in HMDAD, we implemented GRNMFHMDA to predict disease-related microbes and then validated the top 10 predicted microbes by HMDAD or recent literatures.

Asthma, a common long-term inflammatory disease of the airways of the lungs, often starts during childhood and its average number of deaths and death rates (per 100,000 people) respectively reached to 38 and 0.1 in 2016 in the World Health Organization (WHO) European region among 10–14 years old children (Kyu et al., 2018). Here, under the GRNMFHMDA framework, asthma was treated as an investigated disease to explore its potential associated microbes. As a result, 9 out of the top 10 microbes in the prediction list were confirmed to be associated with asthma by experimental literatures (See Table 1). For example, Lactobacillus casei rhamnosus Lcr35, a species of Lactobacillus (1st in the prediction list), was found to be able to attenuate airway inflammation and hyperreactivity in a mouse model of asthma through oral treatment before sensitization (Yu et al., 2010). Besides, Ding et al. (2018) discovered that exosomes derived by Pseudomonas (2nd in the prediction list) aeruginosa could induce protection against allergic sensitization in asthma mice. Another example is that there is a distinct alteration of the sputum microbiota with a greater prominence of Firmicutes (4th in the prediction list) in severe asthma (Zhang et al., 2016).

Obesity, a medical condition in which accumulated excess body fat reaches a certain level that may have a negative effect on health, is a leading preventable cause of death worldwide (Reinier and Chugh, 2015). In recent years, plenty of studies have shown certain associations between obesity and microbes that helps a lot to the prevention and treatment of obesity. For instance, many researchers have demonstrated that methanogens play a specific role in weight gain and the development of obesity in human subjects (Armougom et al., 2009; Krajmalnik-Brown et al., 2012). Not only that, many studies have now been conducted into the potential of probiotics to ameliorate obesity and diabetes (Delzenne et al., 2011; Peterson et al., 2015). Therefore, taking obesity as another investigated disease in the first type of case study, we found that 9 out of the top 10 predicted obesity-related microbes were confirmed by experimental literatures (See Table 2). For the phylum Proteobacteria (1st in the prediction list) which belongs to gram-negative bacteria, the existing study already discovered that it was abundant in the obese group compared with lean group (Park et al., 2015). Besides, as a species of Clostridia (2nd in the prediction list), the presence of Clostridium ramosum in simplified human intestinal (SIHUMI) enhanced diet-induced obesity according to the experiment data of Woting et al. (2014). Moreover, Bacillus, a genus of Clostridia (3rd in the prediction list), was found to have outgrown dramatically in the obesity group by Gao et al. (2018).

More than that, in order to facilitate future researchers to study the disease-related microbes that they are interested in, based on the known associations in HMDAD, we provided the whole prediction list including all pairs between 292 microbes and 39 diseases as well as their predicted association scores (See Supplementary Table 1).

In addition, to prove the predictive applicability of our model on new diseases without known associated microbes, we carried out another case study on a disease via removing all its known associations in HMDAD. In this way, the prediction process of seeking the investigated disease-related microbes can only depend on the information of other known microbe-disease associations (training samples) and the relevant similarity measures. What needs to be emphasized is that only candidate samples (all microbe-disease pairs including the investigated disease) were ranked and then verified in HMDAD. Hence, there was no overlap between training samples and prediction list. In other words, the verification of predicted associations was independent of HMDAD. Type 1 diabetes, a form of diabetes mellitus, is believed to involve a combination of genetic and environmental factors such as dietary agents (Serena et al., 2015), viral infections (Rewers and Ludvigsson, 2016) and gut microbiota (Bibbò et al., 2017). Especially in gut microbiota, the previous study confirmed that the genus Bacteroides is the largest representative of type 1 diabetes-associated dysbiosis that can be modulated by diet (Mejjía-León and Barca, 2015). Thus, considering the significance of studying type 1 diabetes-related microbes, we took type 1 diabetes as the investigated disease to predict its potential associated microbes under the framework of the second type of case study. After implementing GRNMFHMDA, we obtained the ranks of type 1 diabetes' candidate microbes in terms of their association scores (See Table 3). As a result, 8 out of the top 10 predictions were confirmed by HMDAD or recent literatures. For example, Giongo et al. (2011) demonstrated that the Clostridia (1st in the prediction list) sequences increased in control samples (samples of general population) as the abundance of Clostridia decreased overtime in the case samples (samples of patients with type 1 diabetes). Moreover, at the phylum level and at p-values <0.001, Proteobacteria (2nd in the prediction list) was found to be higher in case samples than that in control samples (Brown et al., 2011). Another example is that Lactobacillus strains (a species of Lactobacillus ranking 4th in the prediction list) was found to be able to induce specific changes in the immune system of non-obese diabetic (NOD) mice that can increase or decrease diabetes (Brown et al., 2011).

According to the results presented, GRNMFHMDA consistently achieved an excellent predictive performance in the two types of case studies. With the continuous experimental research on microbe-disease associations, we expect that more and more microbes in the prediction lists generated by our model would be verified in the future.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.