HyMM: hybrid method for disease-gene prediction by integrating multiscale module structure

Abstract

Introduction

The progress of human disease gene discovery has promoted the understanding of the underlying molecular basis of human diseases, but genes known to be associated with diseases only account for a very small proportion of the incidences [1–4]. Traditional approaches such as linkage analysis and genome-wide association studies (GWAS) often provide a long list of candidate genes, requiring expensive and time-consuming experimental identification [5, 6]. Therefore, with the accumulation of biomedical data [7–10], developing computational algorithms for predicting disease-related candidate genes is indispensable to accelerate the discovery of disease-related genes [3, 11, 12].

Organism as a complex biological system is composed of a large number of biomolecules (e.g. genes and proteins) with complex relationships (physical interactions or functional associations), forming a complex biomolecule network system, where the biomolecules exert biological functions through intermolecular synergy while rarely function alone. Human complex diseases can be viewed as the consequences of perturbations or functional abnormalities of associated synergistic biomolecules in the complex network system [13]. Therefore, it is very necessary to study complex diseases and relevant biological phenomena from the perspective of system biology, and biological networks provide an important means for the research of system biology [14–18]. Especially, network-based algorithms have been a popular strategy for the study of disease-related genes [3, 19–28], since genes associated with the same or similar diseases are more similar functionally and their products tend to be highly interconnected in biomolecule networks [4, 29] (see next section). However, how to mine the characteristics of biomolecule networks so as to more effectively explore disease-related genes and related issues is still under continuous exploration due to the inherent complexity of biomolecule networks and the limitation of existing knowledge (e.g. the incompleteness of protein interactome) [15, 17, 30–34].

As we know, module structure as a common property of complex networks is ubiquitous in biomolecule networks, and the modular nature of human diseases can provide useful insights for the study of diseases, but it has not been fully explored in disease-gene prediction [35–37]. Generally, the genes and their products of disease tend to form a disease module due to their high interconnectivity in biomolecule networks [4, 29], but they are usually found to be distributed in multiple modules/subnetworks due to the intrinsic definition of a specific algorithm and the existence of multiscale module structure in the networks [4, 38–41]. The multiscale structure is indeed widespread in biological networks. For example, a module in a protein network may contain several sub-modules, e.g. some protein complexes (such as SAGA) contain several secondary complexes; most of the biological information (e.g. in Gene Ontology) is organized in the form of hierarchical structure. Many algorithms with a flexible resolution parameter have been proposed and applied to mine multiscale module structures in biological networks [42–45], where the resolution parameter can adjust the size or scale of identified modules (see next section). This can provide richer information for studying complex systems such as biomolecule networks, but there are still many challenging issues, such as how to effectively identify the multiscale modules from a network and how to mine the valuable information hidden in the multiscale structure.

To make use of multiscale module structure to more effectively predict disease-related genes, we therefore propose a hybrid method integrating the information of multiscale modules (HyMM) (Figure 1). HyMM extracts a series of module partitions from local to global scales by multiscale modularity optimization (MO) with exponential sampling, and estimates the disease relatedness of genes in partitions by the abundance of disease-related genes in modules. Then, a probabilistic model for integration of gene rankings is designed so as to integrate multiple predictions derived from multiscale module partitions and network propagation, and a parameter estimation strategy based on functional information for Gene Ontology (GO) annotations, pathways or disease genes is proposed to further enhance HyMM’s predictive power.

Figure 1

Open in new tabDownload slide

Workflow of the HyMM method integrating multiscale module structure. (a) Datasets: GG, DD and DG denote the gene–gene, disease–disease and disease–gene associations, respectively; GO and PW denote the GO annotations and pathways of genes, respectively. (b) Extract multiscale module partitions by multiscale algorithm. (c) The extracted module partitions are transformed into matrix representation. (d) A ranking list of genes is generated for each partition matrix; these ranking lists are organized into a ranking matrix of genes, and then, an integrated ranking list of genes based on this ranking matrix is generated by a probabilistic model along with a parameter estimation based on functional information. (e) Generate a ranking list of genes by network propagation. (f) Final rankings of genes are generated by integrating the ranking lists of genes from multiscale modules and network propagation.

The rest of the paper is organized as follows. Firstly, we introduce some related work in this study (including the identification of module structure and the prediction of disease-related genes). Secondly, we present the datasets and details of HyMM, as well as evaluation methods. Thirdly, we study the effectiveness of multiscale module information in disease-gene prediction by combining the functional analysis of multiscale modules; then, by a series of experimental tests, we verify the good performance of HyMM and study the effects of various factors including the definition of conditional probability, multiscale module extraction, parameter estimation based on functional information, sampling of multiscale module partitions and random shuffling of disease-gene associations. Furthermore, we apply HyMM to other datasets as well as specific diseases [e.g. Alzheimer’s disease (AD)] to further demonstrate the effectiveness of HyMM. These results confirm that HyMM can enhance the ability to predict disease-related genes by integrating a multiscale module structure. It is a protocol for disease-gene prediction integrating multiscale module structure, which may become a very useful computational tool for the study of disease-related genes.

Related work

In this study, we focus on the mining of information hidden in the multiscale structure to enhance the ability of disease-gene prediction, which involves two main aspects: disease-gene prediction and (multiscale) module identification (also called community detection or community mining in the field of complex networks [46–48]).

Disease-gene prediction is not only an important issue in computational biology but also is an important field of network medicine/biology [4, 14, 16, 17, 19, 43–45]. Numerous network-based methods for disease-gene prediction have been proposed, based on various approaches, e.g. from homogeneous (HO) network to heterogeneous (HE) network and from single-layer network to multi-layer network [15, 17, 34, 49–52]. For HO network model, for example, Köhler et al. [53] predicted disease-associated genes by the random walk with restart (RWR) on a protein interactome; Chen et al. [54] prioritized disease candidate genes by the k-step markov (KS) method; Hsu et al. [55] developed the gene interconnectedness-based method to rank candidate genes by evaluating the network closeness of them to seeds (known disease-related genes); Zhu et al. [56] proposed the vertex similarity-based (VS) method to discover disease-associated genes. For HE network model, for example, Li et al. [57] proposed the RWR on a disease-gene HE network to infer disease-gene associations; Wu et al. [58] proposed the network-based global inference method called CIPHER to predict human disease-related genes; Xie et al. [59] proposed the bi-random walk (BiRW) to predict disease-gene associations; Singh-Blom et al. [60] predicted disease-gene associations by developing the KATZ measure on a HE network, inspired by social network analyses. For more sophistic network models or techniques, for example, Valdeolivas et al. proposed the RWR on multiplex and HE networks [27]; Xiang et al. [34] proposed the network impulsive dynamics on the multiplex network for disease-gene prediction; Liu et al. [33] proposed a new network embedded representation algorithm to infer pathogenic genes; Xiang et al. [61] proposed a disease-gene-prediction method based on fast network embedding, which can effectively use information in a multi-source HE network constructed by integrating multiple types of association data. Moreover, some module-based algorithms are also applied to the analysis of disease-related genes/modules as well as related issues [42, 62–65]. See references [2, 3, 15, 19, 66–68] for related reviews.

The existing methods in literature have promoted the progress of disease-gene prediction, while the ‘guilt-by-association’ becomes a top-down central principle for predicting disease-related genes in the networks [69]. Evaluating the closeness or distance between candidates and known disease-related genes is a direct strategy to infer disease-related genes in the networks [20, 55, 56], while network propagation (e.g. RWR, KS, RWRH and BiRW) can effectively make use of more information in the whole network to mine potential disease-gene associations [26, 53, 54, 70, 71]. Network propagation (especially the RWR) shows excellent performance in many scenarios [19], so it has been widely applied in the study of bio-entity associations including disease-gene prediction [70].

Biological networks such as protein–protein interaction networks are an important basis for network-based methods in disease-gene prediction, and the mining of biological networks can be helpful for understanding the characteristics of networks and promoting the study of relevant issues. The existence of module structure is an important property of biological networks [72–74], and the research of module structure has been an important topic in the study of complex networks including biological networks. In the past decades, a large number of different types of algorithms based on various approaches (e.g. MO [75], dynamics [76] and statistical inference [77]) were proposed to identify modules (or say communities) in networks, which involve various types of modules (e.g. from single scale to multiple scales, and from non-overlapping to overlapping) and various types of networks (e.g. from unweighted to weighted networks, from undirected to directed networks, and from unsigned to signed networks) [46, 47]. Many of the module identification algorithms, especially MO-based algorithms, have been applied to the study of biological networks, e.g. functional module mining [72, 74, 78], protein complex detection [79–81], and disease module identification [42].

As mentioned above, genes/proteins associated with the same disease tend to form relevant disease modules in a biomolecule network [4, 29], but these genes are usually distributed in multiple modules by specific algorithms [4, 38]. There are several possible reasons for this phenomenon. (i) Complex diseases usually involve functional abnormalities of multiple genes, and these genes may have different functions, playing different roles in the development of complex diseases [82]. (ii) The existing biomolecule networks such as protein–protein interactions are still incomplete [30, 83]. This may cause the detected network modules to be broken and incomplete. (iii) Detected modules in networks are often algorithm-specific, because specific definitions of modules are different for different algorithms [46]. Some algorithms may split a large module into several small submodules in a network, or aggregate several small modules into a large one, because of the existence of a resolution limit that is related to the intrinsic definition or mechanism of algorithms [39–41]. This also implies the existence of a multiscale structure in the network.

In fact, multiscale structure widely exists in various natural and artificial complex networks, including biological networks [84, 85]. In this case, algorithms with flexible resolution parameters, e.g. multiscale MO [86, 87], may more effectively mine the module structure of networks at different scales, where the resolution parameter can be used to tune the scale or size of identified modules [48, 86, 88]. For example, multiscale MO can find relatively large modules in a network when the resolution parameter is small, while it can identify relatively small modules when the resolution parameter is large. This is similar to observing an object from a local to a global scale by a microscope with adjustable resolution parameters. Modules at different scales from local to global ones can be identified by adjusting the resolution parameter in continuous real number space. To study modules at different scales, one generally extracts a set of module partitions corresponding to a set of values sampled from the space of the resolution parameter by a suitable strategy (e.g. exponential sampling).

Multiscale module identification is important for studying biomolecule networks. Dunn et al. [89] have used edge-betweenness clustering to separate protein interaction networks into modules correlating to annotated gene functions, where modules of different sizes can be identified by removing different numbers of edges. Lewis et al. [90] investigated the correlation between the functions of sets of proteins and network module/community structure at multiple resolutions/scales, and they showed that there exist different important scales of module/community structure depending on studied proteins and processes. Wang et al. [91] proposed a fast hierarchical clustering (HC) algorithm using the local metric of edge clustering value, which can uncover the hierarchical organization of functional modules that approximately corresponds to the hierarchical structure of GO annotations. Extended (multiscale) MO was used to identify disease modules [42]. More recently, Zheng et al. developed a multiscale approach called HiDeF to identify robust structures at all scales by integrating the concept of persistent homology with existing community detection algorithms (e.g. multiscale MO).

The mining of multiscale module structure can reveal the features of biological networks at multiple scales, reflecting the correlations of network nodes (e.g. genes) at different levels. This can provide more abundant information for the relevant research involving biological networks, such as network-based disease–gene prediction and protein–protein interaction prediction. Therefore, in this study, we will explore methods for disease-gene prediction by integrating a multiscale module structure.

Materials and methods

Here, we introduce the datasets in this study, and then propose the details of HyMM (including multiscale MO with exponential sampling, disease-relatedness estimation of genes based on multiscale modules, a probabilistic model for integration of multiple gene rankings as well as parameter estimation based on functional information) and evaluation methods. See Figure 1 and Supplementary Note 1 (see Supplementary Data available online at https://academic.oup.com/bib) for the workflow of HyMM.

Datasets

To investigate the predictive ability of algorithms, we employ the disease-gene associations, gene–gene associations and disease-disease associations. In order to conduct functional analysis of modules and parameter estimation based on functional information, we adopt three types of functional groups: GO annotations, PW, and disease-gene sets. See Supplementary Note 2 (see Supplementary Data available online at https://academic.oup.com/bib) for details of datasets.

Disease–gene associations

We use three disease-gene association datasets. (i) The first dataset is an integrated disease-gene dataset [30, 92] retrieved from GWAS and Online Mendelian Inheritance in Man [93]. It is denoted as the Medical Subject Headings Ontology (MeSH) dataset, since MeSH is used to combine the different disease nomenclatures of the two sources into a single standard vocabulary. (ii) The second one is obtained from the DISEASES database, which is a weekly updated web resource for disease-gene associations [94]. (iii) The third one is obtained from the DisGeNet database (https://www.disgenet.org/), which is known as a platform that contains one of the largest publicly available collections of disease-related genes [95]. The UMLS (Unified Medical Language System) diseases in the dataset are mapped into MeSH diseases.

Gene–gene associations

Genes and their products mainly perform their biological functions through their direct or indirect interactions, forming a complex gene–gene association network [83, 96–99]. The gene–gene associations are very important for the study of disease research, since complex diseases are usually considered to be caused by local disturbances of complex biomolecule networks [17, 100, 101]. Here, the gene–gene associations are derived from protein–protein interactions (PPIs). Because single-source protein–protein networks are often incomplete and there exist data noises in existing protein networks, we adopt a comprehensive protein interactome that consists of multiple sources of protein–protein interactions: regulatory interactions, binary interactions from several yeast two-hybrid high-throughput and literature-curated datasets, literature-curated interactions derived mostly from low-throughput experiments, metabolic enzyme-coupled interactions, protein complexes, kinase-substrate pairs and signaling interactions [30]. The network data considers only physical protein interactions with experimental support. The identifiers of proteins are mapped into gene symbols.

These gene–gene associations form a HO network of genes. Furthermore, we construct a disease-gene HE network by integrating gene–gene associations, disease-gene associations and disease-disease associations mentioned above. Note that only the disease-gene associations in the training set are used in the construction of the HE network.

Disease–disease associations

The disease-disease associations can provide useful knowledge for the discovery of disease-related genes. Here, the disease–disease association network is constructed by using the associations between symptoms and diseases. The strengths of these associations between a symptom s$s$ and a disease d$d$ are quantified through the co-occurrence (⁠Cd,s$C_{d,s}$⁠) of their MeSH terms in literature, i.e. the number of PubMed identifiers where two MeSH terms occur together, and then they are normalized as wd,s=Cd,slog(n/ns)$w_{d,s}=C_{d,s}\log (n/n_s)$ by the term frequency-inverse document frequency, where n$n$ denotes the number of diseases and ns≥1$n_s\ge 1$ denotes the number of diseases with symptoms$s$ [102]. Finally, the association score between two diseases is quantified by the cosine similarity scores of their normalized symptom vectors, Score(Vd,Vb)=⟨Vd,Vb⟩/⟨Vd,Vd⟩⟨Vb,Vb⟩−−−−−−−−−−−−√$\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}(V_d,V_b)=⟨V_d,V_b⟩/\sqrt{⟨V_d,V_d⟩⟨V_b,V_b⟩}$⁠, where Vd=(wd,1,wd,2,…,wd,s,…,wd,n)T$V_d={(w_{d,1},w_{d,2},\dots ,w_{d,s},\dots ,w_{d,n})}^T$ denotes the normalized symptom vector of disease d$d$ and ⟨⋅,⋅⟩$⟨\cdot ,\cdot ⟩$ denotes the scalar product of two vectors.

Three types of functional groups

(i) The GO annotations are downloaded from the Molecular Signatures Databases (MSigDB) [103, 104], which omits GO terms with fewer than 5 genes or in very broad categories; (ii) the pathway-gene sets (PW) are also obtained from MSigDB, which were curated from several online pathway databases (such as KEGG and Reactome), publications in PubMed and knowledge of domain experts [105, 106] and (iii) the disease-gene sets (DG) are obtained as mentioned above.

Multiscale MO with exponential sampling

Identification of module/community structure itself is an important issue in the research of networked systems [46, 48, 107]. We here extract module structure from local to global scales by MO with exponential sampling. Moreover, we also consider two other multiscale methods: asymptotic surprise (AS) [88] and fast HC [91]. All of them have flexible resolution parameters to adjust the scale or size of modules. Given a set of reasonably sampled resolution-parameter values, they can generate a set of network module partitions that contain important information of network structure. (see Supplementary Note 3, see Supplementary Data available online at https://academic.oup.com/bib, for details).

Multiscale MO

MO can detect module structure at different scales by varying the resolution parameter γ$\gamma$⁠. It can find global-scale modules (i.e. relatively large modules) when the γ$\gamma$-value is small; it can find local-scale modules (i.e. relatively small modules) when the γ$\gamma$-value is large. It will decompose a network into a set of single-node modules when γ$\gamma$ is large enough, e.g. γ>2M/k2min$\gamma >2M/k_{\mathrm{m}\mathrm{i}\mathrm{n}}^2$⁠, where kmin$k_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is the minimum node degree in the network [109]. This is similar to what happens when we observe objects by a microscope with adjustable resolution. When the resolution of the microscope is high, we can see very local and subtle areas of the object in the field of vision of the microscope; when the resolution is low, we can see its relatively macro and coarse-grained areas, and even its global appearance in the field of vision of the microscope. There is no strict threshold to distinguish between local and global modules, but the limit of ‘local scale’ is that the network is split into a set of single-node modules, and the limit of ‘global scale’ is that the whole network is considered as a large module.

MO needs to be realized with the help of effective optimization algorithms. Here, the Louvain algorithm is applied, because it is a very effective and widely used strategy for optimizing objective functions of module structure in networks, and we have shown that it can be further improved by an effective initialization process and refining process [48, 75, 88].

Exponential sampling of multiscale module partitions

To extract a set of meaningful module partitions from local to global scales {Ψh∣h=1∼H}${\Psi }_h\mid h=1\sim H\}$⁠, where Ψh${\Psi }_h$ denotes the h$h$-th module partition and H$H$ denotes the number of extracted module partitions, we first define a meaningful range of the resolution parameter γ∈[γmin,γmax$\gamma \in [{\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}},{\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$], which covers all possible sampled resolution-parameter values. In a network with N$N$ nodes, γmin${\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}$ and γmax${\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ can theoretically be defined as γmin=max{γ|#[Ψh]=1}${\gamma }_{\mathrm{m}\mathrm{i}\mathrm{n}}=max\{\gamma |\mathrm{\#}[{\Psi }_h]=1\}$ and γmax=min{γ|#[Ψh]=N]}${\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}=min\{\gamma |\mathrm{\#}[{\Psi }_h]=N]\}$⁠, where #[Ψh]$\mathrm{\#}[{\Psi }_h]$ denotes the number of modules in the partition Ψh${\Psi }_h$⁠, but it is usually not easy to obtain accurate interval boundaries. Therefore, we use semi-empirical boundaries according to previous research [87, 110].

According to the set of sampled γ$\gamma$-values, the above multiscale algorithms can extract a set of corresponding module partitions.

Disease-relatedness estimation of genes based on multiscale module structure

Then, the union matrix of gene scorings can be calculated by SG=BSM$S_G=B{S}_M$⁠. Finally, we generate a union matrix of gene rankings, RG=(r→(1)G,r→(2)G,⋯,r→(h)G,⋯,r→(H)G)$R_G=({\overrightarrow{r}}_G^{(1)},{\overrightarrow{r}}_G^{(2)},\cdots ,{\overrightarrow{r}}_G^{(h)},\cdots ,{\overrightarrow{r}}_G^{(H)})$⁠, by the decreasing order of genes’ scores, RG=generank(SG)$R_G=generank(S_G)$⁠. Note that the mean value of ranking is used for genes with the same scores.

According to the above gene scoring/ranking strategy, genes within the same module have the same disease-relatedness scorings/rankings. Therefore, these scoring/ranking lists of genes contain the information of disease relatedness as well as the information of multiscale module structure from low to high resolutions. For example, if a module partition consists of two modules: one contains disease-related genes while not for another, gene scorings/rankings will have two values/levels. If the module with disease-related genes is further split into two sub-modules with disease-related genes, then gene scorings/rankings will have three values/levels if there is no degeneracy of module scorings. The number of values/levels in the scoring/ranking list of genes is closely related to the number of disease-related modules in a module partition. As the resolution increases, we can get the gene scoring/ranking lists with more levels/values, thereby revealing different levels of disease-related information in the network.

This provides another possible way to understand the scorings based on multiscale module partitions.

Probabilistic model for integration of multiple gene rankings

The prior knowledge P(xg)$P(x_g)$ can contribute to the evaluation of the genes’ scores if it is available, or it can be set as a constant if no prior knowledge is available for specific genes.

In order to calculate the final scores of candidate genes, it is necessary to provide an explicit mathematical form of the above conditional probability function (CPF) P(xg|fh)$P(x_g|f_h)$ about each feature (denoted as CPF). For each module partition, we have calculated the scoring list s→(h)G${\overrightarrow{s}}_G^{(h)}$ of candidate genes being related to a disease, and get the ranking list r→(h)${\overrightarrow{r}}^{(h)}$ of the genes. The higher ranking of a gene generally means its larger probability of disease relatedness, which provides a possible way to the definition of CPF. In this study, we design three explicit forms of CPF (denoted as CPF1, CPF2 and CPF3), which correspond to three variants of HyMM.

Without loss of generality, given a ranking list of genes r→(h)G=(r(h)1,r(h)2,…,r(h)g,…,r(h)N)T${\overrightarrow{r}}_G^{(h)}={(r_1^{(h)},r_2^{(h)},\dots ,r_g^{(h)},\dots ,r_N^{(h)})}^T$ based on the decreasing order of disease-relatedness of genes, CPF can be defined as follows,

CPF1: P(xg|fh)∝1−r(h)g/N;$P(x_g|f_h)\propto 1-r_g^{(h)}/N;$

CPF2: P(xg|fh)∝1/r(h)g$P(x_g|f_h)\propto 1/r_g^{(h)}$⁠;

CPF3: P(xg|fh)∝exp(−βhr(h)g)$P(x_g|f_h)\propto \exp (-{\beta }_h{r}_g^{(h)})$⁠,

where β={βh}$\beta =\{{\beta }_h\}$ is a set of tunable parameters (see Supplementary Note 4, see Supplementary Data available online at https://academic.oup.com/bib, for details).

For CPF1, P(xg|fh)$P(x_g|f_h)$ linearly varies with the value of gene ranking; for CPF2, P(xg|fh)$P(x_g|f_h)$ is inversely proportional to the value of gene ranking; for CPF3, P(xg|fh)$P(x_g|f_h)$ exponentially decreases with the value of gene ranking. CPF2 and CPF3 are typically nonlinear functions, which decay more strongly than CPF1 (Figure S1, see Supplementary Data available online at https://academic.oup.com/bib). As a result, CPF2 and CPF3 prefer genes with higher rankings, i.e. genes at the top of the ranking list, because the higher-ranked genes will be assigned the relatively higher possibility values than other lower-ranked genes (Figure S1, see Supplementary Data available online at https://academic.oup.com/bib). This will be conducive to the mining of disease-related genes.

Given the above union matrix of gene rankings from multiscale modules, a comprehensive scoring list s→G${\overrightarrow{s}}_G$ of genes can be generated by the above integration strategy. Then, we get the ranking list of genes r→MM=generank(s→G)${\overrightarrow{r}}_{MM}=generank({\overrightarrow{s}}_G)$ by the decreasing order of genes’ scores (see Figure 1). The above process can be regarded as a raw algorithm for disease-gene prediction (denoted as MM for simplicity).

The scorings/rankings from multiscale modules may provide useful and complementary information for disease-gene prediction, which is different from that of many other algorithms based on various principles, e.g. network propagation. So, we further integrate the ranking list r→MM${\overrightarrow{r}}_{MM}$ with that (denoted by r→TA${\overrightarrow{r}}_{TA}$⁠) of network propagation, e.g. the random walk with a restart in HO/HE networks (see Figure 1). The final scores/rankings of genes will be used to prioritize candidate genes. We will show that this integration can very effectively enhance the ability to predict disease-related genes due to information complementarity.

Parameter estimation based on functional information

Because of the good performance of CPF3, it will be used as the default form of P(xg|fh)$P(x_g|f_h)$ in this study. This integration strategy for multiple gene rankings provides a possible theoretical explanation for classical rank aggregation methods such as Borda count [114], since it degenerates into the arithmetic mean of rankings about multiple features when the parameters β≡1$\beta \equiv 1$⁠. However, different from the classical Borda count, the optimization of β={βh}$\beta =\{{\beta }_h\}$ can further improve the ability to disease-gene prediction, e.g. by using functional information such as GO annotations and PW.

We will use the functional consistency metrics to quantify the functional relevance of network modules and module partitions, based on the set of functional groups for GO/PW/DG. This may provide insights for parameter estimation of β={βh}$\beta =\{{\beta }_h\}$ so as to improve the ability of HyMM to disease-gene prediction.

Evaluation methods

We implement the above procedure of HyMM (including the multiscale algorithms) by Matlab (2016 version). To evaluate the performance of algorithms in disease-gene prediction, we use two evaluation strategies: traditional 5-fold cross-validation (5FCV) and independent test (IndTest). (i) For 5FCV, known disease-related genes for each disease are randomly split into five subsets. In each realization, one of the subsets is treated as a test set, while the rest is treated as a training set. (ii) For IndTest, the disease-gene associations in the MeSH dataset are used as a training set, and the disease-gene associations that only belong to the DisGeNet dataset are used as a test set.

To construct the candidate set of genes, which consists of a test set of genes and a control set of genes, we will construct two kinds of control sets: artificial linkage-interval control set (ALICS) and whole-genome control set (WGCS). (i) For ALICS, each test gene selects 99 control genes from genes closest to this test gene on the same chromosome. This simulates the scenario with disease-related mutation locations (e.g. derived from the genome-wide association study or linkage analysis) [53]. (ii) For WGCS, all unknown genes outside the training and test sets are used as a control set. This simulates the scenario without the information of disease-related mutation locations.

Then, based on the ranking list of candidate genes, several standard evaluation metrics (AUPRC, Recall, and Precision) are used to quantify the performance of prediction algorithms. (i) AUPRC denotes the area under the Precision-Recall curve (PRC), where the PRC curve has Recall on the x-axis and Precision on the y-axis. This is a widely used metric to comprehensively evaluate the performance of algorithms. (ii) Recall measures the ratio of known disease-related genes found in the top-k ranking list compared to the test set, which focuses on how many disease-related genes in the test set have been retrieved. (iii) Precision (Prec) measures the probability of discovering known disease-related genes in the top-k ranking list. Recall and Precision as a function of k-value can provide an intuitive comparison for the local performance of prediction algorithms.

Experimental results

In this section, we first study the effectiveness of multiscale modules and display the good performance of the HyMM framework in disease-gene prediction by a series of experimental tests, including the effect of various factors such as the CPF, multiscale algorithm, parameter estimation based on functional information, and sampling of multiscale module partitions.

Effectiveness of multiscale modules in disease-gene prediction

Here, we study the predictive performance of each module partition being used independently (Figure 2; Figures S2–S4, and Supplementary Note 5, see Supplementary Data available online at https://academic.oup.com/bib). As a whole, the predictive power of the algorithm based on single-scale module partition first has a large upward trend and then a downward trend with the increase of resolution, and the downward trend appears earlier and is more obvious in the HO network. This clearly indicates that different scale module partitions have different levels of importance in disease-gene prediction.

Figure 2

Open in new tabDownload slide

Predictive ability of MO-based module partitions at different scales in disease-gene prediction, as a function of resolution parameter, in HO and HE networks, under different control sets (WGCS and ALICS).

The main reason behind the above phenomenon is that as the resolution increases, modules become more and more fine-grained, and the relevance between nodes in the same module is getting higher and higher. Gradually splitting modules into more fine-grained ones can filter low-relevance nodes while retaining high-relevance nodes in a module, but this will also lead to the loss of information, resulting in the decline of prediction power, because modules without disease-related information are trivial for our scoring strategy.

In order to verify the relationship between the above-mentioned functional relevance of nodes inside modules and the studied scale (resolution), we calculate the functional consistency of module partitions at different scales by using the GO, PW and DG functional groups. The results show that the functional consistency scores of module partitions increase with the increase of resolution (Figure S5, see Supplementary Data available online at https://academic.oup.com/bib). This means that the ratio of similar genes within the modules is increasing with the resolution: these genes are more likely to have the same GO annotations or belong to the same pathway or disease-gene set. This is because the edge density in the modules becomes higher with the increase of resolution. Genes in the modules are more likely to tend to interact with each other and thus have the same or similar functions or participate in a common biological process. Therefore, the module partitions at different scales can provide different levels of information about the relationship between genes. This may provide a more comprehensive understanding of genes and their functions.

Further, with the increase of resolution, disease-related genes also tend to be dispersed into more modules with smaller sizes but stronger functional relevance (e.g. for GO, PW or DG) (Figure S6, see Supplementary Data available online at https://academic.oup.com/bib). As a whole, candidate genes in these modules will be more likely to be disease-related. So, the predictive power of the algorithm based on single-scale module partition gradually increases with the increase of resolution, but the decline of predictive power may appear due to information loss caused by over-filtering of low-relevance nodes, especially when the network is divided into extremely broken modules.

HyMM effectively enhances ability to disease-gene prediction

HyMM outperforms baseline algorithms

Here, we evaluate the performance of HyMM/MM when default setting is used, by comparing to eight baseline algorithms: RWR (Random Walk with Restart) [53], KS (K-Step Markov) [54], VS (Vertex Similarity) [56], ICN (Interconnectedness) [55], RWRH (Random Walk with Restart on Heterogeneous network) [57], CIPHER (Correlating protein Interaction network and PHEnotype network to pRedict disease genes) [58], BiRW (Bi-Random Walk) [59] and KATZ [60] (also see Supplementary Note 6). MM denotes the algorithm that uses only multiscale modules to generate predictive scores. For simplicity, we here define the ratio of performance improvement as (x−y)/y,$(x-y)/y,$where x$x$ and y$y$ denote the results of HyMM and the best baseline algorithm(s), respectively.

The experimental results show that HyMM outperforms all these baseline algorithms, in both the HO and HE networks, under both the ALICS and WGCS control sets (see the results of AUPRC/Recall/Prec in Figure 3; Figures S7 and S10, see Supplementary Data available online at https://academic.oup.com/bib). Specifically, under the ALICS control set, HyMM in the HO network exceeds the best baseline algorithm by 7, 4 and 7% in AUPRC, Recall and Prec metrics, respectively; HyMM in the HE networks exceeds the best baseline algorithm by 27, 25 and 23% in AUPRC, Recall and Prec metrics, respectively. Under the WGCS control set, HyMM in HO exceeds the best baseline algorithms by 9, 21 and 31% in AUPRC, Recall and Prec, respectively; HyMM in HE exceeds the best baseline algorithm by 28, 33 and 32% in AUPRC, Recall and Prec, respectively. The results of top-k Recall/Prec curves have again confirmed the performance of HyMM (Figure 4; Figures S8 and S11, see Supplementary Data available online at https://academic.oup.com/bib).

Figure 3

Open in new tabDownload slide

Performance comparison of HyMM/MM to different baseline algorithms under the ALICS control set. (a and b) AUPRC, (c and d) Recall and (e and f) Precision (Prec) in the HO/HE networks.

Figure 4

Open in new tabDownload slide

Comparison of local performance of HyMM to different baseline algorithms under the ALICS control set. (a and b) Top-k Recall and Precision in the HO network; (c and d) top-k Recall and Precision in the HE network.

HyMM provides a useful framework to enhance ability to disease-gene prediction

We systematically test the performance of the HyMM framework by integrating it with other baseline algorithms. For simplicity, we here define the ratio of performance improvement due to the HyMM framework as (x−y)/y,$(x-y)/y,$where x$x$ and y$y$ denote the results of the improved and original algorithms, respectively. The results show that the HyMM framework can improve the performance of these algorithms in various test scenarios (Figure 5; Figures S9 and S12, see Supplementary Data available online at https://academic.oup.com/bib). For example, under the ALICS control set, the AUPRC, Recall and Prec of CIPHER improve by 87, 88 and 100%, respectively; the AUPRC, Recall and Prec of KATZ improve by 64, 50 and 61%, respectively. Under the WGCS control set, the AUPRC, Recall and Prec of CIPHER improve by 125, 100 and 157%, respectively; the AUPRC, Recall and Prec of KATZ improve by 140, 145 and 158%, respectively. All these results indicate that HyMM is a very effective framework for integrating multiscale modules to promote the ability to disease-gene prediction.

Figure 5

Open in new tabDownload slide

Performance improvement of different baseline algorithms due to the use of multiscale module information in (a) HO and (b) HE networks, under the ALICS control set. −MM and + MM denote the non-use and use of multiscale module information, respectively.

Moreover, ALICS and WGCS simulate the scenarios with and without disease-related mutation locations, respectively. ALICS has a smaller set of candidate genes than WGCS, due to its more information (about mutation locations). Thus, ALICS generally has relatively larger values of evaluation metrics (AUPRC, Recall, and Precision). In fact, this can also be understood from a random point of view, since the probability of randomly selecting correct genes in a smaller candidate set is usually greater.

Comparison of different CPFs

We have compared the performance of HyMM using different CPFs (CPF1/CPF2/CPF3) (Figures S13 and S14, in Supplementary Note 7, see Supplementary Data available online at https://academic.oup.com/bib). In HE, CPF3 can obtain the best performance under both the ALICS and WGCS control sets. In HO, HyMM using CPF3 has comparable or better performance than that using CPF1/CPF2 under the ALICS and WGCS control sets. So, CPF3 is used as the default form of CPF.

Moreover, it is interesting that the nonlinear forms of CPF (e.g. CPF2/CPF3) are better than the linear form (e.g. CPF1). This means that it is beneficial to give more preference to high-ranking genes in this probabilistic framework.

Comparison of different multiscale algorithms

We have compared the performance of HyMM using different multiscale algorithms (MO/AS/HC). Under the ALICS control set, HyMM using MO outperforms HyMM using AS and HC in most cases, while the Recall- and Prec-values of AS in HE are slightly higher than those of MO (Figure 6; Figure S15, in Supplementary Note 8, see Supplementary Data available online at https://academic.oup.com/bib). Under the WGCS control set, HyMM using MO has good performance, although HC is the best in HO, and AS is the best in HE (Figure S16, see Supplementary Data available online at https://academic.oup.com/bib). Moreover, MM using MO is better than that using AS and HC in all the cases. Overall, MO can robustly produce better or comparable performance in various tests, so MO is used as the default choice.

Figure 6

Open in new tabDownload slide

For HyMM/MM, comparison of different multiscale algorithms (MO, AS and HC) under the ALICS control set. (a and b) AUPRC, (c and d) Recall and (e and f) Precision (Prec), in HO and HE networks. HyMM/MM denotes the default algorithms using MO; HyMM-AS/MM-AS and HyMM-HC/MM-HC denote the algorithms using AS and HC, respectively.

Performance improvement through parameter estimation based on functional information

Since module partitions at different scales are of different importance for disease-gene prediction, the optimization of β$\beta$ may further enhance the ability of HyMM, although the default setting has been capable of producing good performance in disease-gene prediction. Due to the close correlation between the functional consistency and predictive power of module partitions, we further use the functional consistency scores of GO/PW/DG as parameter estimation of β$\beta$⁠. The results show that the parameter estimation can indeed improve the ability of HyMM/MM (using MO, AS or HC) comprehensively (Figure 7; Figures S17 and S18, in Supplementary Note 9, see Supplementary Data available online at https://academic.oup.com/bib).

Figure 7

Open in new tabDownload slide

Due to the use of functional information of DG/GO/PW, performance improvement of (a–c) MM and (d–f) HyMM (HE) (using MO, AS and HC), under the ALICS control set. EW denotes that equal weight is used.

Stability to sampling of multiscale module partitions

Sampling of multiscale module partitions is closely related to the number of multiscale module partitions and the amount of information extracted from the network. To study the effect of sampling of multiscale module partitions on prediction performance, we evaluate the performance of HyMM for different values of resolution interval Δlogγ$\Delta log\gamma$⁠. The results show that HyMM using MO/AS is very stable to module partition sampling in various scenarios, while HyMM using HC has relatively large fluctuations and declines (Figure 8; Figures S19 and S20, in Supplementary Note 10, see Supplementary Data available online at https://academic.oup.com/bib).

Figure 8

Open in new tabDownload slide

Performance stability of HyMM using MO/AS/HC (i.e. HyMM_MO, HyMM_AS and HyMM_HC) under the ALICS control set, as a function of resolution interval: (a) AUPRC, (b) Recall and (c) Precision in the HO network; (d) AUPRC, (e) Recall and (f) Precision in the HE network.

Figure 9

Open in new tabDownload slide

Effect of random shuffling of disease-gene associations on predictive performance under the ALICS control set, as a function of ratio of shuffled disease-gene associations: (a and b) AUPRC, (c and d) Recall and (e and f) Precision, in the HO and HE networks, respectively.

Effect of random shuffling of disease-gene associations

Further, we study the effect of random shuffling of disease-gene associations on the predictive performance of algorithms. We generate a series of disease-gene datasets with different degrees of randomization from no randomization to complete randomization by randomly replacing a certain ratio of known disease-gene associations in a dataset with randomly sampled unknown associations; and then we test the performance of algorithms on these datasets (see Supplementary Note 11, see Supplementary Data available online at https://academic.oup.com/bib, for details).

The results show that HyMM consistently outperforms other baseline algorithms with an increasing degree of randomization (Figure 9, Figures S21 and S22, see Supplementary Data available online at https://academic.oup.com/bib). This again confirms the stable and good performance of HyMM in disease-gene prediction. Moreover, as expected, the performance of all algorithms obviously degrades as the degree of randomization increases. This means that HyMM on real datasets is far superior to that on random datasets, and the known disease-gene associations are critical for effectively inferring disease-gene associations. The predictive power of existing prediction algorithms is extremely dependent on the accumulation of confirmed and reliable disease-gene associations, which is the solid basis for the development of disease-gene-prediction algorithms.

Applications to other datasets

In the above sections, we have demonstrated that HyMM has stable and good performance in disease-gene prediction. Here, we further apply HyMM to the other two datasets: the DISEASES and DisGeNET datasets, e.g. by cross-validation and independent test (see Supplementary Note 12, see Supplementary Data available online at https://academic.oup.com/bib).

Performance evaluation in cross-validation

For the DISEASES dataset, the disease terminology in Disease Ontology (DO) database is used, and thus the similarity scores between diseases are calculated based on the DO database by the DOSE package [115]. For the DisGeNet dataset, the UMLS diseases are mapped into the MeSH diseases according to the disease mappings in the DisGeNet database, and thus the MeSH symptom-based disease similarity scores are still used. Here, the disease-gene associations in the two datasets will be used as a benchmark in turn, and we have tested the performance of HyMM in the two datasets by cross-validation experiments (Figures S23 and S24, see Supplementary Data available online at https://academic.oup.com/bib). In the DISEASES dataset, HyMM has higher or comparable values of AUPRC/Recall/Prec than the best baseline algorithm(s) under both the ALICS control set, and HyMM has good overall performance under the WGCS control set (see Supplementary Note 12, see Supplementary Data available online at https://academic.oup.com/bib). In the DisGeNet dataset, HyMM consistently has higher values of AUPRC/Recall/Prec than the best baseline algorithm(s) under both the control sets. These results show that, as in the MeSH dataset, HyMM also has good performance when applied to the datasets, further confirming the effectiveness of HyMM in disease-gene prediction.