A network-based approach to uncover microRNA-mediated disease comorbidities and potential pathobiological implications
The manifestation and clinical severity of human disease are affected by myriad factors, including genetic, epigenetic, lifestyle, and various environmental variables. Identification of disease–disease relationships not only offers insights into disease heterogeneity, but also reveal etiology and pathogenesis of disease comorbidities, , thus driving development of effective therapeutic strategies. , Previous studies designed to map comprehensive disease–disease connections focused mainly on known associations among diseases and associated genes/proteins. However, the predisposition to human disease is dictated by a complex, polygenic, and pleiotropic genetic architecture. Some complex diseases that are mainly driven by environmental or acquired triggers often display more limited genetic risk. Thus, traditional bioinformatics analysis of genetic risk factors offers limited power to detect the true breadth of complex disease–disease relationships. Beyond genetic analysis, shared patterns of gene expression have raised possibilities to inspect disease–disease relationships. Alteration and dysregulation of gene expressions are caused by several biological mechanisms, including microRNA (miRNA) dysregulation. In 1993, Ambros et al. discovered the first type of miRNA (lin-4) in a nematode, revealing for the first time the essential function of miRNA in the posttranscriptional regulation of gene expression. MiRNAs belong to a class of endogenous, small, non-coding RNAs (~22 nucleotides) and play crucial roles in inhibiting the expression of target mRNAs at the posttranscriptional level. Specifically, miRNAs regulate target genes by partially or completely pairing with their 3′ UTR region, thereby reducing the stability of the target miRNA or inhibiting translation to downregulate the expression of genes of interest. This complex regulatory network not only regulates the expression of multiple genes through one miRNA, but also finely regulates the expression of multiple genes by the combination of several miRNAs. Thus, the shared patterns of gene expression regulated by miRNAs may offer possibilities to inspect disease–disease relationships. Currently, more than 30,000 miRNAs within ~200 species have been identified. Cumulative empirical evidences show that miRNAs are closely related to the development, progression, and prognosis of multiple diseases, such as pulmonary vascular disease. , However, it is not obvious whether ascertaining the comprehensive breadth of miRNA-mediated gene networks offer discerning power to reveal important disease–disease relationships. Recent human protein–protein interactome network modeling shows that network-based approaches have raised possibilities to identify disease–disease relationships and drug–disease associations. In this study, we developed a network-based methodology, termed meta-path-based Disease Network (mpDisNet) capturing algorithm, to infer new disease–disease relationships from miRNA-mediated network perspectives. We built a heterogeneous miRNA–gene–disease network by assembling four biological networks: disease–miRNA, miRNA–gene, gene–disease, and the human protein–protein interactome (Table ). Specifically, mpDisNet searches a specific meta-path (a meta-path is a path linking two specified nodes in a network mode) based on a Random Walk algorithm to reconstruct the heterogeneous neighbors of a node. Specifically, we utilized a heterogeneous skip-gram model to solve the network representation of the nodes in mpDisNet (Fig. ). We found that mpDisNet displayed a higher performance in inferring disease–disease relationships compared with traditional miRNA-overlapping approaches. Via t-distributed stochastic neighbor embedding (t-SNE) analysis, the reduced dimension graphs generated by the disease–miRNA–gene and disease–gene networks reveal that mpDisNet can effectively distinguish different class of human diseases, offering potential pathobiological implications. We further identified pulmonary disease comorbidities (e.g., lung cancer-asthma and asthma-chronic obstructive pulmonary disease) with potential miRNA-mediated pathobiological mechanisms. If broadly applied, mpDisNet would offer a powerful network-based tool for identification of disease–disease relationships for multiple complex diseases from heterogeneous biological networks. 

Pipeline of mpDisNet MpDisNet infers miRNA-mediated disease–disease relationships based on the topology of multiple networks among diseases, miRNAs, and genes (Fig. ). The pipeline of mpDisNet has four key steps (see Methods section): (i) network data integration: we reconstructed a heterogenous network by assembling four experimentally validated networks, including disease–miRNA, miRNA–gene, disease–gene, and the human interactome networks (Table ); (ii) meta-path-based Random Walks: we reconstructed heterogeneous neighbors of the nodes using the random walk of the meta-path and generated instance sequences; (iii) heterogeneous skip-gram: we generated the multidimensional vector for each disease by the skip-gram from the instance sequences; and (iv) network-based inferring disease–disease relationships: we calculated the disease–disease cosine similarities based on the multidimensional vectors generated from the skip-gram (iii). The detailed pipeline of mpDisNet is illustrated in Fig. . Performance of mpDisNet We compared mpDisNet with miRNA-overlap measure on the experimentally validated disease–miRNA association network (see Methods section). Herein, mpDisNet is the result of selecting the meta-path M1 (disease–miRNA–gene–gene–miRNA–disease) and M3 (disease–gene–gene–disease) in an integrated heterogeneous network (Fig. ). For miRNA-overlap measure, we assume that the set of miRNAs corresponding to disease A is A m , and the corresponding set of disease B is B m . We calculated disease–disease similarity based on overlap measure as below: 1 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$S_{{}} = }  B_{}}}{{A_{}  B_{}}}$$ S overlap = A m ⋂ B m A m ⋃ B m We selected the top 300 pairs of the highest similarity disease pairs (Supplementary Table ) obtained by miRNA-overlap measure and mpDisNet, and plotted two network graphs of miRNA-overlap measure (Fig. ) and mpDisNet (Fig. ), respectively. The node color of each disease is classified according to the disease pathobiological classification from a previous study. Overall, the mpDisNet (Fig. ) can capture clinically reported disease–disease comorbidities in the same pathobiological categories of specific diseases, outperforming miRNA-overlap measure (Fig. ). For example, associations among obesity (Mesh ID: D009765), diabetes mellitus (Mesh ID: D003920), cystic fibrosis (Mesh ID: D003550), osteoporosis (Mesh ID: D010024), and metabolic syndrome X (Mesh ID: D024821) are well captured by mpDisNet (Fig. ). For cardiovascular disease, the significant associations among heart disease (myocardial infarction), coronary artery disease, atherosclerosis, ischemia, and hypertension are successfully identified by mpDisNet as well (Fig. ). For neurological diseases, the mpDisNet-predicted relationships among schizophrenia, bipolar disorder, and Alzheimer’s disease were consistent with a recent study. Finally, multiple types of cancer are found to share a strong association identified by mpDisNet, consistent with recent pan-cancer studies. , Altogether, mpDisNet identifies potentially well-known disease–disease relationships. To validate performance of mpDisNet further, we collected 220 clinically reported disease–disease pairs from a previous study. We found that these 220 disease–disease pairs can be correctly re-identified by mpDisNet. However, miRNA-overlap measure can only identify 120 pairs. We plotted the network map (Fig. ) of mpDisNet-predicted 100 comorbid disease pairs (Supplementary Table ) which are not identified by miRNA-overlap measure. For example, mpDisNet successfully identifies the associations of autoimmune lymphoproliferative syndrome with bipolar disorder, cataract, celiac disease, and Crohn disease. In addition, cerebral infarction is associated with several diseases or syndromes, including friedreich ataxia, long QT Syndrome, multiple endocrine neoplasia Type 1, osteogenesis imperfecta, retinitis pigmentosa, telangiectasia, hereditary hemorrhagic, and thalassemia, identified by mpDisNet as well (Fig. and Supplementary Table ). We next turned to evaluate the receiver operating characteristic (ROC) and precision-recall curves based on 66 clinically reported disease–disease pairs (Supplementary Table ) derived from the previously published implicit semantic similarity measure. We found that mpDisNet showed a reasonable accuracy (the area under ROC [AUROC = 0.65] and the area under precision-recall curve [AUPR] = 0.68, Fig. ) in inferring the clinically reported disease–disease pairs, outperforming that of miRNA-overlap measure (AUROC = 0.59 and AUPR = 0.56, Fig. ). In addition, mpDisNet showed a reasonable accuracy (AUROC = 0.67 and AUPR = 0.66) in inferring the clinically reported disease–disease pairs on an external validation set, revealing high generalizability. Altogether, mpDisNet reveals high accuracy in inferring disease–disease relationships, outperforming traditional miRNA-overlap measure. Biological interpretation of mpDisNet We next turned to investigate whether the underlying miRNA-mediated subnetworks identified by mpDisNet can offer potential pathobiological mechanisms for the inferred disease–disease relationships. Specifically, we integrated two networks into a single heterogeneous network and evaluated two meta-paths M1 (disease–miRNA–gene–gene–miRNA–disease) and M3 (disease–gene–gene–disease) as shown in Fig. . The multidimensional vectors of the two meta-paths were obtained by random walk and skip-gram, and then the multidimensional vectors were concatenated to infer disease–disease relationships (see Methods). We then performed dimensionality reduction visualization analysis using a t-SNE algorithm. We removed diseases with unknown classification and kept diseases with well-known pathobiological annotations with at least seven types of diseases in each category. In the dimensionality reduction diagram (Fig. ), a closer distance between two diseases reveals a higher relevant pathobiological relationship. We found that the same pathobiological categories of diseases are clustered by the multidimensional vectors (Fig. ), indicating that the underlying miRNA-mediated pathobiological pathways can be identified by mpDisNet. Network-based identification of miRNA-mediated pathobiological pathways between lung cancer and asthma As shown in Fig. , we found a strong association of cancers (e.g., lung neoplasms) with asthma and COPD. This finding is consistent with recent meta-analyses, suggesting the potential associations of COPD and asthma with several cancer types such as lung cancer. , For example, shortness of breath and respiratory distress often increase the suffering of advanced-stage lung cancer patients. , However, the underlying disease pathways for lung cancer-associated asthma remain unclear. Asthma is a condition characterized by chronic inflammation of the lungs, including airway hyper-reactivity, excessive mucous formation, and respiratory obstruction. We asserted that lung cancer-associated asthma may be caused from tumor cell microenvironments, such as cross-talk pro-inflammatory pathway. For example, recent studies showed that micro-environmental inflammation by tumor cell-immune cell cross-talk may induce lung cancer-associated pulmonary hypertension. , We therefore performed a multi-layer human interactome network analysis to inspect the miRNA-mediated pathobiological pathways for lung cancer-associated asthma via mpDisNet (Fig. ). For example, two highlighted miRNAs, hsa-mir-7a and hsa-mir-155, play important roles in both lung cancer , and asthma, , which are involved in multiple meta-paths in Fig. . Hsa-mir-34a was reported as a tumor suppressor gene by inhibiting non-small cell lung cancer (NSCLC) growth and suppressing the CD44hi stem-like NSCLC cells. , We found that a meta-path of hsa-mir-34a-SAA1-APBB1 may involve in the lung cancer-associated asthma by meta-path-based network analysis within the human protein–protein interactome (Fig. ). SSA1 , encoding serum amyloid A1, activates the NLRP3 inflammasome and promotes asthma in mice. Thus, hsa-mir-34a that mediates lung tumor growths, may involve in inflammasome-mediated pathways in asthma as well. We next examined whether we can identify novel miRNA-mediated pathways for lung cancer-associated asthma. Figure reveals that a meta-path of hsa-mir-17-STK11/LKB1 plays a key role in lung cancer by regulating cancer cell metabolism. – STK11/LKB1 is a central regulator of T cell development, activation and metabolism. In addition, the T cell plays an important functional role in asthma as well. Collectively, hsa-mir-17-STK11/LKB1 may offer a potential pathobiological pathway for lung cancer-associated asthma. In summary, potential miRNA-mediated disease pathways captured by mpDisNet offer candidate biomarkers in understanding of pathobiological mechanisms of lung cancer-associated asthma. However, these candidate network biomarkers identified by mpDisNet are warranted by experimental or clinical validation further. Network-based identification of miRNA-mediated pathobiological pathways between COPD and asthma Asthma and COPD are obstructive pulmonary diseases that have affected millions of people all over the world. They are two diseases with differences in etiology, symptoms, type of airway inflammation, inflammatory cells, mediators, consequences of inflammation, response to therapy and course. The similarities in airway inflammation in severe asthma and COPD and good response to combination therapies in both diseases suggest that they may share some pathophysiologic characteristics. , We next turned to inspect the miRNA-mediated pathways between asthma-COPD. Both hsa-let-7a (differentially expressed in patients with severe asthma ) and hsa-let-7b play important roles in asthma by targeting pro-inflammatory pathways. We found two meta-paths, including hsa-let-7a-CASP3-CCND1-hsa-mir-20a and hsa-let-7b-CCND2-FOXO4-hsa-mir-499a between asthma and COPD, via mpDisNet (Fig. ). Genetic studies and in vitro observations have shown potential associations of CCND1 and CCND2 with asthma and COPD. – In addition, CASP3 was reported to play a functional role in airway epithelial apoptosis , and pro-inflammatory cytokines (FOXO4) may contribute to regulation of muscle atrophy and smooth muscle cell migration. , Altogether, miRNA-mediated airway epithelial apoptosis and pro-inflammatory cytokine pathways (hsa-let-7a and hsa-let-7b) may offer potential mechanisms for the overlapping syndrome between asthma and COPD. In addition, several mpDisNet-predicted meta-paths, such as hsa-mir-148b-ADAM33-PGD-hsa-mir-1 and hsa-mir-221-ACTB-BUB1-hsa-mir-196a (Fig. ) may offer new pathobiological pathways to explain the asthma-COPD comorbidity as well. – 

MpDisNet infers miRNA-mediated disease–disease relationships based on the topology of multiple networks among diseases, miRNAs, and genes (Fig. ). The pipeline of mpDisNet has four key steps (see Methods section): (i) network data integration: we reconstructed a heterogenous network by assembling four experimentally validated networks, including disease–miRNA, miRNA–gene, disease–gene, and the human interactome networks (Table ); (ii) meta-path-based Random Walks: we reconstructed heterogeneous neighbors of the nodes using the random walk of the meta-path and generated instance sequences; (iii) heterogeneous skip-gram: we generated the multidimensional vector for each disease by the skip-gram from the instance sequences; and (iv) network-based inferring disease–disease relationships: we calculated the disease–disease cosine similarities based on the multidimensional vectors generated from the skip-gram (iii). The detailed pipeline of mpDisNet is illustrated in Fig. . 

We compared mpDisNet with miRNA-overlap measure on the experimentally validated disease–miRNA association network (see Methods section). Herein, mpDisNet is the result of selecting the meta-path M1 (disease–miRNA–gene–gene–miRNA–disease) and M3 (disease–gene–gene–disease) in an integrated heterogeneous network (Fig. ). For miRNA-overlap measure, we assume that the set of miRNAs corresponding to disease A is A m , and the corresponding set of disease B is B m . We calculated disease–disease similarity based on overlap measure as below: 1 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$S_{{}} = }  B_{}}}{{A_{}  B_{}}}$$ S overlap = A m ⋂ B m A m ⋃ B m We selected the top 300 pairs of the highest similarity disease pairs (Supplementary Table ) obtained by miRNA-overlap measure and mpDisNet, and plotted two network graphs of miRNA-overlap measure (Fig. ) and mpDisNet (Fig. ), respectively. The node color of each disease is classified according to the disease pathobiological classification from a previous study. Overall, the mpDisNet (Fig. ) can capture clinically reported disease–disease comorbidities in the same pathobiological categories of specific diseases, outperforming miRNA-overlap measure (Fig. ). For example, associations among obesity (Mesh ID: D009765), diabetes mellitus (Mesh ID: D003920), cystic fibrosis (Mesh ID: D003550), osteoporosis (Mesh ID: D010024), and metabolic syndrome X (Mesh ID: D024821) are well captured by mpDisNet (Fig. ). For cardiovascular disease, the significant associations among heart disease (myocardial infarction), coronary artery disease, atherosclerosis, ischemia, and hypertension are successfully identified by mpDisNet as well (Fig. ). For neurological diseases, the mpDisNet-predicted relationships among schizophrenia, bipolar disorder, and Alzheimer’s disease were consistent with a recent study. Finally, multiple types of cancer are found to share a strong association identified by mpDisNet, consistent with recent pan-cancer studies. , Altogether, mpDisNet identifies potentially well-known disease–disease relationships. To validate performance of mpDisNet further, we collected 220 clinically reported disease–disease pairs from a previous study. We found that these 220 disease–disease pairs can be correctly re-identified by mpDisNet. However, miRNA-overlap measure can only identify 120 pairs. We plotted the network map (Fig. ) of mpDisNet-predicted 100 comorbid disease pairs (Supplementary Table ) which are not identified by miRNA-overlap measure. For example, mpDisNet successfully identifies the associations of autoimmune lymphoproliferative syndrome with bipolar disorder, cataract, celiac disease, and Crohn disease. In addition, cerebral infarction is associated with several diseases or syndromes, including friedreich ataxia, long QT Syndrome, multiple endocrine neoplasia Type 1, osteogenesis imperfecta, retinitis pigmentosa, telangiectasia, hereditary hemorrhagic, and thalassemia, identified by mpDisNet as well (Fig. and Supplementary Table ). We next turned to evaluate the receiver operating characteristic (ROC) and precision-recall curves based on 66 clinically reported disease–disease pairs (Supplementary Table ) derived from the previously published implicit semantic similarity measure. We found that mpDisNet showed a reasonable accuracy (the area under ROC [AUROC = 0.65] and the area under precision-recall curve [AUPR] = 0.68, Fig. ) in inferring the clinically reported disease–disease pairs, outperforming that of miRNA-overlap measure (AUROC = 0.59 and AUPR = 0.56, Fig. ). In addition, mpDisNet showed a reasonable accuracy (AUROC = 0.67 and AUPR = 0.66) in inferring the clinically reported disease–disease pairs on an external validation set, revealing high generalizability. Altogether, mpDisNet reveals high accuracy in inferring disease–disease relationships, outperforming traditional miRNA-overlap measure. 

We next turned to investigate whether the underlying miRNA-mediated subnetworks identified by mpDisNet can offer potential pathobiological mechanisms for the inferred disease–disease relationships. Specifically, we integrated two networks into a single heterogeneous network and evaluated two meta-paths M1 (disease–miRNA–gene–gene–miRNA–disease) and M3 (disease–gene–gene–disease) as shown in Fig. . The multidimensional vectors of the two meta-paths were obtained by random walk and skip-gram, and then the multidimensional vectors were concatenated to infer disease–disease relationships (see Methods). We then performed dimensionality reduction visualization analysis using a t-SNE algorithm. We removed diseases with unknown classification and kept diseases with well-known pathobiological annotations with at least seven types of diseases in each category. In the dimensionality reduction diagram (Fig. ), a closer distance between two diseases reveals a higher relevant pathobiological relationship. We found that the same pathobiological categories of diseases are clustered by the multidimensional vectors (Fig. ), indicating that the underlying miRNA-mediated pathobiological pathways can be identified by mpDisNet. 

As shown in Fig. , we found a strong association of cancers (e.g., lung neoplasms) with asthma and COPD. This finding is consistent with recent meta-analyses, suggesting the potential associations of COPD and asthma with several cancer types such as lung cancer. , For example, shortness of breath and respiratory distress often increase the suffering of advanced-stage lung cancer patients. , However, the underlying disease pathways for lung cancer-associated asthma remain unclear. Asthma is a condition characterized by chronic inflammation of the lungs, including airway hyper-reactivity, excessive mucous formation, and respiratory obstruction. We asserted that lung cancer-associated asthma may be caused from tumor cell microenvironments, such as cross-talk pro-inflammatory pathway. For example, recent studies showed that micro-environmental inflammation by tumor cell-immune cell cross-talk may induce lung cancer-associated pulmonary hypertension. , We therefore performed a multi-layer human interactome network analysis to inspect the miRNA-mediated pathobiological pathways for lung cancer-associated asthma via mpDisNet (Fig. ). For example, two highlighted miRNAs, hsa-mir-7a and hsa-mir-155, play important roles in both lung cancer , and asthma, , which are involved in multiple meta-paths in Fig. . Hsa-mir-34a was reported as a tumor suppressor gene by inhibiting non-small cell lung cancer (NSCLC) growth and suppressing the CD44hi stem-like NSCLC cells. , We found that a meta-path of hsa-mir-34a-SAA1-APBB1 may involve in the lung cancer-associated asthma by meta-path-based network analysis within the human protein–protein interactome (Fig. ). SSA1 , encoding serum amyloid A1, activates the NLRP3 inflammasome and promotes asthma in mice. Thus, hsa-mir-34a that mediates lung tumor growths, may involve in inflammasome-mediated pathways in asthma as well. We next examined whether we can identify novel miRNA-mediated pathways for lung cancer-associated asthma. Figure reveals that a meta-path of hsa-mir-17-STK11/LKB1 plays a key role in lung cancer by regulating cancer cell metabolism. – STK11/LKB1 is a central regulator of T cell development, activation and metabolism. In addition, the T cell plays an important functional role in asthma as well. Collectively, hsa-mir-17-STK11/LKB1 may offer a potential pathobiological pathway for lung cancer-associated asthma. In summary, potential miRNA-mediated disease pathways captured by mpDisNet offer candidate biomarkers in understanding of pathobiological mechanisms of lung cancer-associated asthma. However, these candidate network biomarkers identified by mpDisNet are warranted by experimental or clinical validation further. 

Asthma and COPD are obstructive pulmonary diseases that have affected millions of people all over the world. They are two diseases with differences in etiology, symptoms, type of airway inflammation, inflammatory cells, mediators, consequences of inflammation, response to therapy and course. The similarities in airway inflammation in severe asthma and COPD and good response to combination therapies in both diseases suggest that they may share some pathophysiologic characteristics. , We next turned to inspect the miRNA-mediated pathways between asthma-COPD. Both hsa-let-7a (differentially expressed in patients with severe asthma ) and hsa-let-7b play important roles in asthma by targeting pro-inflammatory pathways. We found two meta-paths, including hsa-let-7a-CASP3-CCND1-hsa-mir-20a and hsa-let-7b-CCND2-FOXO4-hsa-mir-499a between asthma and COPD, via mpDisNet (Fig. ). Genetic studies and in vitro observations have shown potential associations of CCND1 and CCND2 with asthma and COPD. – In addition, CASP3 was reported to play a functional role in airway epithelial apoptosis , and pro-inflammatory cytokines (FOXO4) may contribute to regulation of muscle atrophy and smooth muscle cell migration. , Altogether, miRNA-mediated airway epithelial apoptosis and pro-inflammatory cytokine pathways (hsa-let-7a and hsa-let-7b) may offer potential mechanisms for the overlapping syndrome between asthma and COPD. In addition, several mpDisNet-predicted meta-paths, such as hsa-mir-148b-ADAM33-PGD-hsa-mir-1 and hsa-mir-221-ACTB-BUB1-hsa-mir-196a (Fig. ) may offer new pathobiological pathways to explain the asthma-COPD comorbidity as well. – 

Understanding of disease–disease relationships is important for the diagnosis, prevention, and treatment of the human disease. Most of the existing comorbid data are from the medical records analysis of clinical patients. This method requires a large amount of data calculation and has many interference factors. Recent remarkable development of systems biology technologies and network medicine approaches raised possibilities to predict disease comorbidities from human protein–protein interactome. , In order to integrate biological networks to predict disease–disease relationships, we presented a network-based methodology, termed mpDisNet, to infer disease–disease relationships from miRNA regulatory network perspective. Specifically, we constructed a comprehensive, multi-layer biological network connecting diseases, miRNA, and genes. We employed a skip-gram algorithm to obtain the multidimensional feature vectors of disease and then calculated the disease–disease similarities from the reduced informative multidimensional vectors. We demonstrated that mpDisNet can identify both clinically reported and new disease–disease associations, outperforming miRNA-overlap measure. Moreover, mpDisNet offers miRNA-mediated pathobiological pathways by searching miRNA meta-paths from the human protein–protein interactome, as we showcased for lung cancer-associated asthma and asthma-COPD. However, comprehensive validation for more mpDisNet-predicted disease–disease relationships are warranted in the future. We highlighted several significant contributions in the current study. We assembled four comprehensive networks, including disease–miRNA, miRNA–genes, disease–gene, and the human protein–protein interactome to search the meta-paths by mpDisNet. In this way, we can utilize the complementary information from different biological networks compared with traditional network-based approaches using single type of data. , Network analysis further shows that integrating miRNA-mediated network can improve the capability in inferring disease–disease relationships, offering a new network-based tool for assessment of disease comorbidities. In addition, the network-based framework presented in mpDisNet could be applied for prediction of drug–target interactions, gene–gene (protein–protein) interactions, RNA–RNA interactions, and other biological networks as well. Finally, the new disease–disease relationships inferred by mpDisNet may offer potential candidate network biomarkers for better understanding of underlying pathobiological pathways from miRNA network perspective. We acknowledged several potential limitations in current network-based framework of mpDisNet. First, when the known miRNA associated with disease is fewer, the comorbidity between disease pairs computed by miRNA-mediated networks may be false positive. Second, potential literature data bias (e.g., degree/connectivity of well-studied miRNAs/proteins) may generate a potential false positive rate. Third, each random walk requires a specific meta-path, and the choice of this single meta-path may also affect performance of mpDisNet. In the future, we may improve mpDisNet by integrating more comprehensive biological networks, analyzing the relevant associations in tissue-specific networks in which the disease occurs, adopting more flexible random walk strategies. In summary, this study offers a network-based, systems biology methodology for comprehensive identification of disease–disease relationships from miRNA regulatory network perspective. From a translational perspective, if broadly applied, mpDisNet would offer a powerful network-based tool for understanding of clinical comorbidities for multiple complex diseases from heterogeneous biological networks, a significant challenge of precision medicine. 

Reconstruction of heterogeneous networks We reconstructed a heterogenous miRNA–gene–disease network by assembling four types of networks: (a) disease–miRNA, (b) miRNA–gene, (c) disease–gene, and (d) the human protein–protein interactome networks. Disease–miRNA network We collected experimentally validated disease–miRNA associations from two databases: miR2Disease and HMDD v3.0. All disease terms were annotated by Medical Subject Headings (MeSH) and Unified Medical Language System (UMLS) vocabularies. The disease–miRNA associations in two databases were combined and the duplicate associations were removed. Finally, we kept a total of 7669 associations connecting 691 miRNAs with 394 diseases in this study. miRNA–gene network We collected the known miRNA targets to build miRNA–gene networks from miRTarBase database. We annotated all protein-coding genes using gene Entrez ID, chromosomal location, and the official gene symbols from the National Center for Biotechnology Information (NCBI) database. In this study, we only kept the data from Homo sapiens. After excluding duplicate associations, 163,090 miRNA–gene associations connecting 568 miRNAs with 14,762 human genes were used. Disease–gene network We assembled disease–gene associations from four public databases: the Online Mendelian Inheritance in Man (OMIM), HuGE Navigator, PharmGKB, and Comparative Toxicogenomics Database (CTD). All disease terms were annotated using MeSH vocabularies, and the genes were annotated using the Entrez IDs and official gene symbols from the NCBI database. Duplicated pairs from different data sources were deleted. In total, we obtained 50,589 disease–gene associations connecting 2684 genes with 394 unique disease terms. The human protein–protein interactome To build a comprehensive human protein–protein interactome, we focused on high-quality protein–protein interactions (PPIs) with five types of experimental evidences: (i) Binary PPIs tested by high-throughput yeast-two-hybrid (Y2H) systems; , (ii) Kinase-substrate interactions by literature-derived low-throughput and high-throughput experiments; (iii) Literature-curated PPIs identified by affinity purification followed by mass spectrometry (AP-MS), Y2H and by literature-derived low-throughput experiments; (iv) PPIs from protein three-dimensional (3D) structures; and (v) Signaling networks supported by literature-derived low-throughput experiments. The genes were mapped to their Entrez ID based on the NCBI database as well as their official gene symbols based on GeneCards ( http://www.genecards.org/ ). Duplicated PPIs and all computationally predicted data, such as evolutionary analysis, metabolic associations, and gene co-expression data, were deleted. The resulting updated human interactome used in this study includes 246,995 PPIs connecting 16,706 unique proteins. The detailed descriptions are provided in our recent studies. , Meta-path-based random walks We employed a meta-path-based random walk to capture the semantic and structural correlation between different types of nodes. Given a heterogeneous network, G = ( V , E , F ), and meta-path, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$P:V_1 {  }^{R_1} V_2 {  }^{R_2} V_3 {  }^{R_3}  V_f {  }^{R_f} V_{f + 1}   {  }^{R_{l - 1}} V_l$$ P : V 1 → R 1 V 2 → R 2 V 3 → R 3 ⋯ V f → R f V f + 1 ⋯ → R l − 1 V l , the transition probability in step i was defined as follows: 2 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$P({v^{i + 1}}  v_f^i,p) = \{ {}{}(v_f^i)}}}&{({v^{i + 1}},v_f^i)  E, ({v^{i + 1}}) = f + 1} \\0&{({v^{i + 1}},v_f^i)  E, ({v^{i + 1}})\,  \,f + 1} \\0&{({v^{i + 1}},v_f^i) \;\; E}} .$$ P ( v i + 1 ∣ v f i , p ) = 1 ∣ N f + 1 ( v f i ) ∣ ( v i + 1 , v f i ) ∈ E , ∅ ( v i + 1 ) = f + 1 0 ( v i + 1 , v f i ) ∈ E , ∅ ( v i + 1 ) ≠ f + 1 0 ( v i + 1 , v f i ) ∉ E where [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$v_f^i  V_f$$ v f i ∈ V f , and [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$N_{f + 1}(v_f^i)$$ N f + 1 ( v f i ) represent the set of nodes belonging to the type, V f +1 , in the neighborhood of node, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$v_f^i$$ v f i . In other words, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$v^{i + 1}  V_{f + 1}$$ v i + 1 ∈ V f + 1 , walking is on the condition of a preset meta-path, P . Moreover, meta-paths are generally used on symmetric paths, that is, its first node type V 1 is the same with the last one V l , facilitating its recursive for random walks, i.e., 3 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$P( {v^{i + 1}{}v_f^i} ) = p( {v^{i + 1}{}v_l^i} ),\,{}\,f = l$$ P v i + 1 ∣ v f i = p v i + 1 ∣ v l i , if f = l The meta-path-based random walk strategy ensures that the semantic relationships among different types of nodes are properly conserved in the reconstructed heterogeneous network. Heterogeneous skip-gram Furthermore, we employed a heterogeneous skip-gram representation learning model. The heterogeneous skip-gram is a modification based on the original Skip-gram model, by adding the superposition of different node types. For a heterogeneous network, G = ( V , E , F ), each node, ν , and each edge, e , are associated with their mapping functions, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$ ( v ):V  F_V(| {F_V} |\, > \,1)$$ φ v : V → F V ( F V > 1 ) and [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$ ( e ):E  F_E$$ ψ e : E → F E , respectively. Given a node, ν , maximizes the probability that the heterogeneous context, N f ( ν ), [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$( {f  F_V} )$$ f ∈ F V is as follows: 4 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${}_{}  { }_{v V}  { }_{f F_V}  { }_{c_f N_f(v)} {}\,p( {c_f{}v; } )$$ argmax θ ∑ v ϵ V ∑ f ϵ F V ∑ c f ϵ N f ( v ) log p c f ∣ v ; θ where N f ( ν ) denotes the neighborhood of ν with the f th type of nodes. The conditional probability, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$p( {c_f{}v; } )$$ p c f ∣ v ; θ , is defined as a softmax function and adjusted to a specific node type, f , as follows: 5 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$p( {c_f{}v; } ) =  X_v}}}{{ { }_{u_f V_f} e^{X_{u_f} X_v}}}$$ p c f ∣ v ; θ = e X c f ⋅ X v ∑ u f ϵ V f e X u f ⋅ X v where X v is the v th row of X , which is the embedding vector for node v ; V f represents the node type set of type, f , in the network. This specifies a multinomial distribution for each type in the output layer of the last layer of skip-gram. According to the negative sampling in Word2vec, the above function is defined as follows: 6 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${O}(X) =   ( {X_{c_f}  X_v} ) +  {}_{m = 1}^M { _{u_f^m P_f( {u_f} )}} [ {{} ( { - X_{u_f^m}  X_v} )} ]$$ O ( X ) = log σ X c f ⋅ X v + ∑ m = 1 M E u f m ~ P f u f log σ − X u f m ⋅ X v where [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${}( x ) = }}$$ σ x = 1 1 + e − x and P f ( u f ) are pre-defined distributions by the type of node of neighbor, c f , that aims to predict from which a negative node [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$u_f^m$$ u f m is drawn from for M times. The gradients of the above pre-defined distributions are derived as follows: 7 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$O(X)}{{ X_{u_f^m}}} = ( { ( {X_{u_f^m}  X_v - {}_{c_f}[u_f^m]} )} )X_v$$ ∂ O ( X ) ∂ X u f m = σ X u f m ⋅ X v − I c f [ u f m ] X v 8 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$(X)}}{{ X_v}} =  { }_{m = 0}^M ( { ( {X_{u_f^m}  X_v - { {I}}_{c_f}[u_f^m]} )} )X_{u_f^m}$$ ∂ O ( X ) ∂ X v = ∑ m = 0 M σ X u f m ⋅ X v − I c f [ u f m ] X u f m where [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${ I}_{c_f}[u_f^m]$$ I c f [ u f m ] is an indicator function to indicate whether [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$u_f^m$$ u f m is the neighborhood context node c f . When m = 0, then [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$u_f^0 = c_f$$ u f 0 = c f . The model is optimized by using the stochastic gradient descent algorithm. Network-based inferring disease–disease relationships The network-based similarities between two diseases can be calculated based on single meta-path or multiple meta-paths. In this study, we evaluated three meta-paths (M1, M2, M3) to infer disease–disease relationships. For M1 (disease–miRNA–gene–gene–miRNA–disease) as shown in Fig. , we randomly walked in disease–miRNA–gene heterogeneous network based on meta-path M1 for 50 steps. Each walk includes 251 nodes. We run 1000 random walks for each disease and 1000 random walk instance sequences are generated. By inputting all the sequences into heterogeneous skip-gram, we obtained the representation vectors of each disease. Then, we calculated the cosine similarity between diseases based on these vectors. In this way, we calculated the disease similarity for meta-path M2 (disease–miRNA–gene–gene–gene–miRNA–disease), M3 (disease–gene–gene–disease) as well. We predicted disease–disease relationships based on multiple meta-paths by concatenating the representation vectors learned from each meta-path and then calculated the cosine similarity between the concatenated vectors. Therefore, we assembled a disease–miRNA–gene network and a disease–gene network into a heterogeneous network. In this integrated heterogeneous network, we selected the meta-paths M1 and M3, respectively. The multidimensional vectors of the two meta-paths can be obtained by random walk and skip-gram, and then the multidimensional vectors were concatenated to infer disease–disease relationships. The detailed network-based analyses are provided in our recent studies. , , 

We reconstructed a heterogenous miRNA–gene–disease network by assembling four types of networks: (a) disease–miRNA, (b) miRNA–gene, (c) disease–gene, and (d) the human protein–protein interactome networks. Disease–miRNA network We collected experimentally validated disease–miRNA associations from two databases: miR2Disease and HMDD v3.0. All disease terms were annotated by Medical Subject Headings (MeSH) and Unified Medical Language System (UMLS) vocabularies. The disease–miRNA associations in two databases were combined and the duplicate associations were removed. Finally, we kept a total of 7669 associations connecting 691 miRNAs with 394 diseases in this study. miRNA–gene network We collected the known miRNA targets to build miRNA–gene networks from miRTarBase database. We annotated all protein-coding genes using gene Entrez ID, chromosomal location, and the official gene symbols from the National Center for Biotechnology Information (NCBI) database. In this study, we only kept the data from Homo sapiens. After excluding duplicate associations, 163,090 miRNA–gene associations connecting 568 miRNAs with 14,762 human genes were used. Disease–gene network We assembled disease–gene associations from four public databases: the Online Mendelian Inheritance in Man (OMIM), HuGE Navigator, PharmGKB, and Comparative Toxicogenomics Database (CTD). All disease terms were annotated using MeSH vocabularies, and the genes were annotated using the Entrez IDs and official gene symbols from the NCBI database. Duplicated pairs from different data sources were deleted. In total, we obtained 50,589 disease–gene associations connecting 2684 genes with 394 unique disease terms. The human protein–protein interactome To build a comprehensive human protein–protein interactome, we focused on high-quality protein–protein interactions (PPIs) with five types of experimental evidences: (i) Binary PPIs tested by high-throughput yeast-two-hybrid (Y2H) systems; , (ii) Kinase-substrate interactions by literature-derived low-throughput and high-throughput experiments; (iii) Literature-curated PPIs identified by affinity purification followed by mass spectrometry (AP-MS), Y2H and by literature-derived low-throughput experiments; (iv) PPIs from protein three-dimensional (3D) structures; and (v) Signaling networks supported by literature-derived low-throughput experiments. The genes were mapped to their Entrez ID based on the NCBI database as well as their official gene symbols based on GeneCards ( http://www.genecards.org/ ). Duplicated PPIs and all computationally predicted data, such as evolutionary analysis, metabolic associations, and gene co-expression data, were deleted. The resulting updated human interactome used in this study includes 246,995 PPIs connecting 16,706 unique proteins. The detailed descriptions are provided in our recent studies. , 

We collected experimentally validated disease–miRNA associations from two databases: miR2Disease and HMDD v3.0. All disease terms were annotated by Medical Subject Headings (MeSH) and Unified Medical Language System (UMLS) vocabularies. The disease–miRNA associations in two databases were combined and the duplicate associations were removed. Finally, we kept a total of 7669 associations connecting 691 miRNAs with 394 diseases in this study. 

We collected the known miRNA targets to build miRNA–gene networks from miRTarBase database. We annotated all protein-coding genes using gene Entrez ID, chromosomal location, and the official gene symbols from the National Center for Biotechnology Information (NCBI) database. In this study, we only kept the data from Homo sapiens. After excluding duplicate associations, 163,090 miRNA–gene associations connecting 568 miRNAs with 14,762 human genes were used. 

We assembled disease–gene associations from four public databases: the Online Mendelian Inheritance in Man (OMIM), HuGE Navigator, PharmGKB, and Comparative Toxicogenomics Database (CTD). All disease terms were annotated using MeSH vocabularies, and the genes were annotated using the Entrez IDs and official gene symbols from the NCBI database. Duplicated pairs from different data sources were deleted. In total, we obtained 50,589 disease–gene associations connecting 2684 genes with 394 unique disease terms. 

To build a comprehensive human protein–protein interactome, we focused on high-quality protein–protein interactions (PPIs) with five types of experimental evidences: (i) Binary PPIs tested by high-throughput yeast-two-hybrid (Y2H) systems; , (ii) Kinase-substrate interactions by literature-derived low-throughput and high-throughput experiments; (iii) Literature-curated PPIs identified by affinity purification followed by mass spectrometry (AP-MS), Y2H and by literature-derived low-throughput experiments; (iv) PPIs from protein three-dimensional (3D) structures; and (v) Signaling networks supported by literature-derived low-throughput experiments. The genes were mapped to their Entrez ID based on the NCBI database as well as their official gene symbols based on GeneCards ( http://www.genecards.org/ ). Duplicated PPIs and all computationally predicted data, such as evolutionary analysis, metabolic associations, and gene co-expression data, were deleted. The resulting updated human interactome used in this study includes 246,995 PPIs connecting 16,706 unique proteins. The detailed descriptions are provided in our recent studies. , 

We employed a meta-path-based random walk to capture the semantic and structural correlation between different types of nodes. Given a heterogeneous network, G = ( V , E , F ), and meta-path, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$P:V_1 {  }^{R_1} V_2 {  }^{R_2} V_3 {  }^{R_3}  V_f {  }^{R_f} V_{f + 1}   {  }^{R_{l - 1}} V_l$$ P : V 1 → R 1 V 2 → R 2 V 3 → R 3 ⋯ V f → R f V f + 1 ⋯ → R l − 1 V l , the transition probability in step i was defined as follows: 2 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$P({v^{i + 1}}  v_f^i,p) = \{ {}{}(v_f^i)}}}&{({v^{i + 1}},v_f^i)  E, ({v^{i + 1}}) = f + 1} \\0&{({v^{i + 1}},v_f^i)  E, ({v^{i + 1}})\,  \,f + 1} \\0&{({v^{i + 1}},v_f^i) \;\; E}} .$$ P ( v i + 1 ∣ v f i , p ) = 1 ∣ N f + 1 ( v f i ) ∣ ( v i + 1 , v f i ) ∈ E , ∅ ( v i + 1 ) = f + 1 0 ( v i + 1 , v f i ) ∈ E , ∅ ( v i + 1 ) ≠ f + 1 0 ( v i + 1 , v f i ) ∉ E where [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$v_f^i  V_f$$ v f i ∈ V f , and [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$N_{f + 1}(v_f^i)$$ N f + 1 ( v f i ) represent the set of nodes belonging to the type, V f +1 , in the neighborhood of node, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$v_f^i$$ v f i . In other words, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$v^{i + 1}  V_{f + 1}$$ v i + 1 ∈ V f + 1 , walking is on the condition of a preset meta-path, P . Moreover, meta-paths are generally used on symmetric paths, that is, its first node type V 1 is the same with the last one V l , facilitating its recursive for random walks, i.e., 3 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$P( {v^{i + 1}{}v_f^i} ) = p( {v^{i + 1}{}v_l^i} ),\,{}\,f = l$$ P v i + 1 ∣ v f i = p v i + 1 ∣ v l i , if f = l The meta-path-based random walk strategy ensures that the semantic relationships among different types of nodes are properly conserved in the reconstructed heterogeneous network. 

Furthermore, we employed a heterogeneous skip-gram representation learning model. The heterogeneous skip-gram is a modification based on the original Skip-gram model, by adding the superposition of different node types. For a heterogeneous network, G = ( V , E , F ), each node, ν , and each edge, e , are associated with their mapping functions, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$ ( v ):V  F_V(| {F_V} |\, > \,1)$$ φ v : V → F V ( F V > 1 ) and [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$ ( e ):E  F_E$$ ψ e : E → F E , respectively. Given a node, ν , maximizes the probability that the heterogeneous context, N f ( ν ), [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$( {f  F_V} )$$ f ∈ F V is as follows: 4 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${}_{}  { }_{v V}  { }_{f F_V}  { }_{c_f N_f(v)} {}\,p( {c_f{}v; } )$$ argmax θ ∑ v ϵ V ∑ f ϵ F V ∑ c f ϵ N f ( v ) log p c f ∣ v ; θ where N f ( ν ) denotes the neighborhood of ν with the f th type of nodes. The conditional probability, [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$p( {c_f{}v; } )$$ p c f ∣ v ; θ , is defined as a softmax function and adjusted to a specific node type, f , as follows: 5 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$p( {c_f{}v; } ) =  X_v}}}{{ { }_{u_f V_f} e^{X_{u_f} X_v}}}$$ p c f ∣ v ; θ = e X c f ⋅ X v ∑ u f ϵ V f e X u f ⋅ X v where X v is the v th row of X , which is the embedding vector for node v ; V f represents the node type set of type, f , in the network. This specifies a multinomial distribution for each type in the output layer of the last layer of skip-gram. According to the negative sampling in Word2vec, the above function is defined as follows: 6 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${O}(X) =   ( {X_{c_f}  X_v} ) +  {}_{m = 1}^M { _{u_f^m P_f( {u_f} )}} [ {{} ( { - X_{u_f^m}  X_v} )} ]$$ O ( X ) = log σ X c f ⋅ X v + ∑ m = 1 M E u f m ~ P f u f log σ − X u f m ⋅ X v where [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${}( x ) = }}$$ σ x = 1 1 + e − x and P f ( u f ) are pre-defined distributions by the type of node of neighbor, c f , that aims to predict from which a negative node [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$u_f^m$$ u f m is drawn from for M times. The gradients of the above pre-defined distributions are derived as follows: 7 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$O(X)}{{ X_{u_f^m}}} = ( { ( {X_{u_f^m}  X_v - {}_{c_f}[u_f^m]} )} )X_v$$ ∂ O ( X ) ∂ X u f m = σ X u f m ⋅ X v − I c f [ u f m ] X v 8 [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$(X)}}{{ X_v}} =  { }_{m = 0}^M ( { ( {X_{u_f^m}  X_v - { {I}}_{c_f}[u_f^m]} )} )X_{u_f^m}$$ ∂ O ( X ) ∂ X v = ∑ m = 0 M σ X u f m ⋅ X v − I c f [ u f m ] X u f m where [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$${ I}_{c_f}[u_f^m]$$ I c f [ u f m ] is an indicator function to indicate whether [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$u_f^m$$ u f m is the neighborhood context node c f . When m = 0, then [12pt]{minimal}
				
				 
				 
				 
				
				
				
				
				$$u_f^0 = c_f$$ u f 0 = c f . The model is optimized by using the stochastic gradient descent algorithm. 

The network-based similarities between two diseases can be calculated based on single meta-path or multiple meta-paths. In this study, we evaluated three meta-paths (M1, M2, M3) to infer disease–disease relationships. For M1 (disease–miRNA–gene–gene–miRNA–disease) as shown in Fig. , we randomly walked in disease–miRNA–gene heterogeneous network based on meta-path M1 for 50 steps. Each walk includes 251 nodes. We run 1000 random walks for each disease and 1000 random walk instance sequences are generated. By inputting all the sequences into heterogeneous skip-gram, we obtained the representation vectors of each disease. Then, we calculated the cosine similarity between diseases based on these vectors. In this way, we calculated the disease similarity for meta-path M2 (disease–miRNA–gene–gene–gene–miRNA–disease), M3 (disease–gene–gene–disease) as well. We predicted disease–disease relationships based on multiple meta-paths by concatenating the representation vectors learned from each meta-path and then calculated the cosine similarity between the concatenated vectors. Therefore, we assembled a disease–miRNA–gene network and a disease–gene network into a heterogeneous network. In this integrated heterogeneous network, we selected the meta-paths M1 and M3, respectively. The multidimensional vectors of the two meta-paths can be obtained by random walk and skip-gram, and then the multidimensional vectors were concatenated to infer disease–disease relationships. The detailed network-based analyses are provided in our recent studies. , , 

Supplementary Tables 1-4 nr-reporting-summary 

