The first decade of Genome Wide Association Studies (GWAS) has uncovered a wealth of disease-associated variants. Two important derivations will be the translation of this information into a multiscale understanding of pathogenic variants and leveraging existing data to increase the power of existing and future studies through prioritization. We explore edge prediction on heterogeneous networks—graphs with multiple node and edge types—for accomplishing both tasks. First we constructed a network with 18 node types—genes, diseases, tissues, pathophysiologies, and 14 MSigDB (molecular signatures database) collections—and 19 edge types from high-throughput publicly-available resources. From this network composed of 40,343 nodes and 1,608,168 edges, we extracted features that describe the topology between specific genes and diseases. Next, we trained a model from GWAS associations and predicted the probability of association between each protein-coding gene and each of 29 well-studied complex diseases. The model, which achieved 132-fold enrichment in precision at 10% recall, outperformed any individual domain, highlighting the benefit of integrative approaches. We identified pleiotropy, transcriptional signatures of perturbations, pathways, and protein interactions as influential mechanisms explaining pathogenesis. Our method successfully predicted the results (with AUROC = 0. 79) from a withheld multiple sclerosis (MS) GWAS despite starting with only 13 previously associated genes. Finally, we combined our network predictions with statistical evidence of association to propose four novel MS genes, three of which (JAK2, REL, RUNX3) validated on the masked GWAS. Furthermore, our predictions provide biological support highlighting REL as the causal gene within its gene-rich locus. Users can browse all predictions online (http: //het. io). Heterogeneous network edge prediction effectively prioritized genetic associations and provides a powerful new approach for data integration across multiple domains.
In the last decade, genome-wide association studies (GWAS) have been established as the main strategy to map genetic susceptibility in dozens of complex diseases and phenotypes. Despite the success of this approach in mapping variation in thousands of loci to hundreds of complex phenotypes [1–5], researchers are now confronted with the challenge of maximizing the scientific contribution of existing GWAS datasets, whose undertakings represented a substantial investment of human and monetary resources from the community at large [6]. A central assumption in GWAS is that every region in the genome (and hence every gene) is a-priori equally likely to be associated with the phenotype in question. As a result, small effect sizes and multiple comparisons limit the pace of discovery. However, rational prioritization approaches may afford an increase in study power while avoiding the constraints and expense related to expanded sampling. One such way forward is the current trend of analyzing the combined contribution of susceptibility variants in the context of biological pathways, rather than single SNPs [7]. For example, Yaspan et al described an approach that aggregates variants of interest from a GWAS into biological pathways using genomic randomization to control for multiple testing and minimize type I error [8]. The popular software PLINK also includes an option to evaluate groups of associations at the gene level, thus enabling pathway analysis by computing enriched gene sets [9]. A less explored but potentially revealing strategy is the integration of diverse sources of data to build more accurate and comprehensive models of disease susceptibility. Several strategies have been attempted to identify the mechanisms underlying pathogenesis and use these insights to prioritize genes for genetic association analyses. Gene-set enrichment analyses identify prevalent biological functions amongst genes contained in disease-associated loci [10,11]. Gene network approaches search for neighborhoods of genes where disease-associated loci aggregate [12,13]. Jia et al. reported dmGWAS, a strategy to integrate association signals from GWAS into the human protein interaction network [14]. A similar approach was developed by our group and tested in two large studies comprising more than 15,000 cases [15]. Literature mining techniques aim to chronicle the relatedness of genes to identify a subset of highly-related associated genes. For example, Raychaudhuri et al. reported the Gene Relationships Among Implicated Loci (GRAIL) algorithm, an approach to assess relationships among genomic disease regions by text mining of PubMed abstracts [16]. Prioritization strategies generally rely on user-provided loci as the sole input and do not incorporate broader disease-specific knowledge. Typically, the proportion of genome-wide significant discoveries in a given GWAS is low, thus leaving little high-confidence signal for seed-based approaches to build from. To overcome this limitation, here we aimed at characterizing the ability of various information domains to identify pathogenic variants across the entire compendium of complex disease associations. Using this multiscale approach, we developed a framework to prioritize both existing and future GWAS analyses and highlight candidate genes for further analysis. To approach this problem, we resorted to a method that integrated diverse information domains naturally. Heterogeneous (or multipartite) networks are a class of networks which contain multiple types of entities (nodes) and relationships (edges or links), and provide a data structure capable of expressing diversity in an intuitive and scalable fashion. Most existing techniques available for network analysis have been developed for homogeneous networks [17–19] and are not directly extensible to heterogenous networks. Accordingly, the early network analyses in disease biology concentrated on homogeneous networks. However in the last half-decade, the complexity of biological systems has spurred interest in heterogeneous approaches. While still a developing field, network-based biological data integration has been pursued using a variety of techniques. Approaches such as GeneMANIA, weight and then project individual data sources onto a single dimension, enabling homogeneous network algorithms to be used to characterize the resulting graphs [20–22]. Other techniques operate on multi-relational (single node type, multiple edge type) networks, for example by taking into account relationships among local clusters and considering the full topology of weighted gene association networks [13,20,23]. Bipartite networks contain two node types and therefore work well for predicting relationships between entities of two different types (such as disease-gene associations or drug-disease indications) following a ‘guilt-by-association’ paradigm [24–26]. Other approaches incorporate greater-dimension heterogeneous networks as input but conflate types and, while improving predictions compared to simpler approaches, cannot effectively identify influential network components [27,28]. Heterogeneous networks of arbitrary complexity have also been applied for edge prediction without a formalized feature extraction methodology, which requires manual descriptor determination for each new network design [29]. Recently, new types of edge prediction methods were reported that naturally accommodate any size heterogeneous network. These include data fusion by matrix-factorization [30–33] and metapath-based techniques [34,35]. This type of intermediate data fusion can treat all data sources directly (i. e. without transforming data into “disease space”) and has been successfully used to infer disease similarities [31] and predict gene function in slime mold and baker’s yeast [32]. A metapath-based approach was recently developed by researchers studying social sciences to predict future coauthorship [34] and provides an intuitive framework and interpretable models and results. An advantage of metapath-based approaches is that they preserve the network structure and provide the flexibility to explore a diverse set of descriptors. In this work, we extended this methodology to predict the probability that an association between a gene and disease exists.
Using publicly-available databases and standardized vocabularies, we constructed a heterogeneous network with 40,343 nodes and 1,608,168 edges (Fig 1, S2 Data). Databases were selected based on quality, reusability, throughput, and their aggregate ability to provide a diversified, multiscale portrayal of biology (see Resource selection in Methods). The network was designed to encode entities and relationships relevant to pathogenesis. The network contained 18 node types (metanodes) and 19 edge types (metaedges), displayed in Fig 2A. Entities represented by metanodes consisted of diseases, genes, tissues, pathophysiologies, and gene sets for 14 MSigDB collections [36,37] including pathways [38,39], perturbation signatures, motifs [40,41], and Gene Ontology (GO) domains [42] (Table 1). Relationships represented by metaedges consisted of gene-disease association, disease pathophysiology, disease localization, tissue-specific gene expression, protein interaction, and gene-set membership for each MSigDB collection (Table 2). Gene-disease associations were extracted from the GWAS Catalog [5] by overlapping associations into disease-specific loci. Loci were classified as low or high-confidence based on p-value and sample size of the corresponding GWAS (see Associations in Methods and S1 Fig). When possible, for each loci, the most-commonly reported gene across studies was designated as primary and subsequently considered responsible for the association. Additional genes reported for the loci were considered secondary. Only high-confidence primary associations were included in the network yielding 938 associations between 99 diseases and 711 genes (S2 Fig visualizes a subset of these associations). To describe the network topology connecting a specific gene and disease, we computed 24 features, each describing a different aspect of connectivity. Each feature corresponds to a type of path (metapath [53]) originating in a given source gene and terminating in a given target disease. The biological interpretation of a feature derives from its metapath (S1 Table), and features simply quantify the prevalence of a specific metapath between any gene-disease pair. To quantify metapath prevalence, we adapted an existing method originally developed for social network analysis (PathPredict) [34], and developed a new metric called degree-weighted path count (DWPC, Fig 2D), which we employed in all but two features. The DWPC downweights paths through high-degree nodes when computing metapath prevalence. The strength of downweighting depends on a single parameter (w), which we optimized to w = 0. 4 and that outperformed the top metric resulting from PathPredict (S3 Fig) [34]. We calculated DWPC features for the 22 metapaths of length 3 or less that originated with a gene and terminated with disease. Two non-DWPC features were included to assess the pleiotropy of the source gene and the polygenicity of the target disease. Referred to as ‘path count’ features, they respectively equal the number of diseases associated with the source gene and the number of genes associated with the target disease. For all features, paths with duplicate nodes were excluded, and, if present, the association edge between the source gene and target disease was masked. Further analysis focused on the 29 diseases with at least ten associated genes (Table 3). The 698 high-confidence primary associations of these 29 diseases were considered positives—gene-disease pairs with positive experimental relationships (as defined in the Associations section of Methods, S2 Fig). The remaining 551,823 (i. e. unassociated) gene-disease pairs were considered negatives. Low-confidence or secondary associations were excluded from either set. We partitioned gene-disease pairs into training (75%) and testing (25%) sets and created a training network with the testing associations removed. To learn the importance of each feature and model the probability of association of a given gene-disease pair, we used regularized logistic regression which is designed to prevent overfitting and accurately estimate regression coefficients when models include many features. Elastic net regression is a regression method that balances two regularization techniques: ridge (which performs coefficient shrinkage) and lasso (which performs coefficient shrinkage and variable selection) [54]. On the training set, we optimized the elastic net mixing parameter, a single parameter behind the DWPC metric, and two edge-inclusion thresholds (S3 Fig). While cross-validated performance was similar across elastic net mixing parameters, ridge demonstrated the greatest consistency (S3 Fig), and thus we proceeded with logistic ridge regression as the primary model for predictions. We extracted network-based features for gene-disease pairs from the training network and modeled the training set. We next evaluated performance on the 25% of gene-disease pairs (175 positives, 137,956 negatives) withheld for testing. Our predictions achieved an area under the ROC curve (AUROC) of 0. 83 (Fig 3A) demonstrating an excellent performance in retrieving hidden associations. Importantly, we did not observe any significant degradation of performance from training to testing (Fig 3A), indicating that our disciplined regularization approach avoided overfitting and that predictions for associations included in the network were not biased by their presence in the network. Furthermore, we observed that at 10% recall (the classification threshold where 10% of true positives were predicted as positives), our predictions achieved 16. 7% precision (the proportion of predicted positives that were correct). Since the prevalence of positives in our dataset was 0. 13%, the observed precision represents a 132-fold enrichment over the expected probability under a uniform distribution of priors (as in GWAS). As a next step in our analysis, we recomputed features on the complete network, which now included the previously withheld testing associations. On all positives and negatives, we fit a ridge model (the primary model for predictions) and a lasso model (for comparison). Standardized coefficients (Fig 4) indicate the effect attributed to each feature by the models. The lasso highlighted features that captured pleiotropy (4 features), pathways (2), transcriptional signatures of perturbations (1) and protein interactions (1). Despite the parsimony of the lasso, performance was similar between models with training AUROCs of 0. 83 (ridge) and 0. 82 (lasso). However, since multiple features from a correlated group may be causal, the lasso model risks oversimplifying. Ridge regression disperses an effect across a correlated group of features, providing users greater flexibility when interpreting predictions. From the ridge model, we predicted the probability that each protein-coding gene was associated with each analyzed disease (S1 Data) and built a webapp to display the predictions (http: //het. io/disease-genes/browse). Using Markov chain randomized edge-swaps, we created 5 permuted networks. Since metaedge-specific node degree is preserved, features extracted from the permuted network retain unspecific effects. These effects include general measures a disease’s polygenicity and a gene’s pleiotropy, multifunctionality, and tissue-specificity. On the first permuted network, we partitioned associations into training and testing sets. Testing associations were masked from the network, features were computed, and a ridge model was fit on the training gene-disease pairs. Compared to the unpermuted-network model, testing performance was noticeably inferior: the AUROC declined from 0. 83 (Fig 3A) to 0. 79 and the AUPRC (area under the precision-recall curve) declined from 0. 06 (Fig 3B) to 0. 02 (S4 Fig). We interpret the modest decline in AUROC but marked reduction in AUPRC as a direct consequence of the permutation’s particularly detrimental effect on top predictions (S4 Fig). In other words, edge-specificity was crucial for top predictions, while general effects gleaned from node degree performed reasonably well when ranking the entire spectrum of protein-coding genes for association. A commonly-overlooked finding is that the discriminatory ability of gene networks largely relies on node-degree rather than the edge-specificity [55]. However, we found that for top predictions—which are the only predictions considered by many applications—edge-specificity was critical. Interestingly, predictions from the permuted-network model displayed a reduced dynamic range with none exceeding 4%, while predictions from the unpermuted-network model exceeded 75% (S4 Fig). Therefore, even though they achieve reasonable AUROC, the permuted-network predictions would have little utility as prior probabilities in a bayesian analysis where dynamic range is crucial. Furthermore, the signal present in permuted-network features was greatly diminished: few features survived the lasso’s selection resulting in an average lasso AUROC of 0. 70 versus 0. 80 for ridge (S5 Fig). Permuting the network significantly reduced the predictiveness of features based on pleiotropy (2 features), protein interactions (2), transcriptional signatures of perturbations (1), tissue-specificity (1), pathways (3), and immunologic signatures (1) (S2 Table). Six of the eight features selected by the lasso and eight of the top ten ridge features (ranked by standardized coefficients) were negatively affected by the permutation. Since our modeling technique preferentially selected/weighted features affected by permutation, we can infer that network components where edge-specificity matters underlie a large portion of predictions. We assessed the informativeness of each feature by calculating feature-specific AUROCs. Feature-specific AUROCs universally exceeded 0. 5, indicating that network connectivity, regardless of type, positively discriminates associations. However, performance varied widely by feature and within feature from disease to disease (Fig 5). Top performing domains consisted of transcriptional signatures of perturbations (AUROC = 0. 74), immunologic signatures (0. 70), and pleiotropy (0. 68,0. 67,0. 64,0. 63). Notably, the models greatly outperformed any individual feature, highlighting the importance of an integrative approach. Features whose metapaths originate with an association (GaD) metaedge measure pleiotropy (S1 Table). The four pleiotropic features were among the top performing features that did not rely on set-based gene categorization (Fig 5). Of the four features, GaD (any disease) had the highest AUROC despite its lack of disease-specificity, reflecting both the sparsity of disease-specific features and the existence of genetic overlap between seemingly disparate diseases. GaDmPmD and GaDaGaD performed best for immunologic diseases and were affected by permutation, indicating that genetic overlap was greatest between immunologic diseases. On the other hand, the performance of GaDlTlD did not decrease after permutation indicating disease colocalization was not a primary driver of genetic overlap (S5 Fig and S2 Table). We also observed that the lasso regression model discarded the majority of features with a minimal performance deficit, suggesting redundancy among features. Indeed, pairwise feature correlations showed moderate collinearity among features (S6 Fig). Collinearity was especially pervasive with respect to the Perturbations feature, explaining its threefold increase in standardized coefficient in the lasso versus ridge model. The disappearance of all but one other MSigDB-based feature in the lasso model indicated that Perturbations—the feature traversing chemical and genetic transcriptional signatures of perturbations—exhausted meaningful gene-set characterization. In other words, the faulty molecular processes behind pathogenesis align with and are encapsulated by the processes perturbed by chemical and genetic modifications. The Immunologic signatures feature—traversing gene-sets characterizing “cell types, states, and perturbations within the immune system”—was highly predictive and correlated with Perturbations. As expected this feature performed best for diseases with an immune pathophysiology. The one well-performing neoplastic disease (Fig 5) was chronic lymphocytic leukemia, a hematologic cancer with a strong immune component [56]. Additionally, the performance of both the Perturbation and Immunologic features was affected by permutation indicating information beyond the extent of a gene’s multifunctionality was encoded. Existing network-based gene-prioritization methods, frequently rely solely on protein-protein interactions. Our results supported incorporating protein interactions as the two interactome-based features were discriminatory (AUROCs = 0. 65,0. 56) and affected by permutation. However, when compared to the integrative models or other top-performing features, performance of features that relied solely on the interactome was severely limited. Pathways, another founding resource for many approaches, proved important with KEGG selected by the lasso and all three pathway resources (AUROCs = 0. 61 for KEGG, 0. 60 for Reactome, 0. 55 for BioCarta) affected by permutation. The GeTlD feature—measuring to what extent a gene is expressed in tissues affected by the disease in question—peaked in performance around AUROC = 0. 58 (S3 Fig), was affected by permutation, and required no preexisting knowledge of associated genes. In other words, while approaches based on tissue-specificity may have limited predictive ability on their own, they are broadly applicable (i. e. less susceptible to knowledge bias) and provide orthogonal information that could enhance the overall performance of a model. For each type of gene set, we evaluated the effect of increased sparsity on performance by randomly subsampling gene set nodes or edges and measuring the resulting AUROC of the affected feature (S7 Fig). Robustness refers to a gene set collection’s ability to withstand a high extent of masking with little performance deficit. Several of the top gene sets had this property, especially GO processes (where supersets are common), which may indicate nodal redundancy. Contrastingly, the MSigDB gene set with the fewest nodes, KEGG, experienced a more immediate and linear decline in performance. Since KEGG avoids duplication and is stringently and manually curated, this finding is expected. To investigate whether the high predictivity of certain gene set collections was due only to size, we compared performance when subsampling nodes to the KEGG level (crosses in S7 Fig). The two top performing collections, perturbations and immunologic signatures, which also happen to be large, continued to perform better than the majority of complete collections. While performance benefited from increasing densities, a resource’s sparsity often reflects an intrinsic property of the underlying information type. Therefore, when identifying influential mechanisms of pathogenesis, we prefered unadjusted comparisons using the complete network. The WTCCC2 multiple sclerosis (MS) GWAS tested 465,434 SNPs for 9,772 cases and 17,376 controls and identified over 50 independently associated loci [57]. Since the GWAS Catalog excludes targeted arrays (such as ImmunoChip), this study remains the largest MS GWAS in the Catalog. To evaluate our method’s ability to prioritize associations identified in a future study, we masked the WTCCC2 MS study from the GWAS Catalog and created a pre-WTCCC2 network. The number of high-confidence primary MS associations was thus reduced from 50 to 13, with the 37 novel genes identified by WTCCC2 available to evaluate performance. On the pre-WTCCC2 network, we extracted features, fit a ridge model, and predicted each gene’s probability of association with MS. Amongst all 18,993 potentially novel genes, the 37 WTCCC2 genes were ranked highly (AUROC = 0. 79, Fig 6). Finally, we designed a framework for discovering and validating novel MS genes that incorporates our network-based predictions. Meta2. 5 is a meta-analysis of all MS GWAS prior to the WTCCC2 study [58]. We calculated genewise p-values for Meta2. 5 using VEGAS [59] and observed a large enrichment in nominally significant (p < 0. 05) genes, suggesting multiple potential associations (S9 Fig). We combined this set of experimental candidates with the top predictions from the pre-WTCCC2 network to discover genes with both strong statistical and biological evidence of association (S12 Data). To ensure novelty, we excluded genes from GWAS-established MS loci and the extended MHC region. We chose a threshold (S3 Table) for network-based predictions that performed well in prioritizing the genes identified by WTCCC2 (Fig 6). This strategy discovered four genes, three of which—JAK2, REL, RUNX3—achieved Bonferroni validation on VEGAS-converted WTCCC2 p-values (Table 4). The probability of the observed validation rate occurring under random prioritization is 0. 01 (S3 Table), demonstrating that incorporating our network-based predictions as a prior increased study power. JAK2 displays overexpression in MS-affected Th17 cells [60] and was implicated in an interactome-based prioritization of GWAS [15]. RUNX3, a transcription factor influencing T lymphocyte development, has been associated with celiac disease [61] and ankylosing spondylitis [62] and was hypermethylated in systemic lupus erythematosus patients [63]. The region containing REL was uncovered in a recent MS ImmunoChip-based study with 14,498 cases [64]. For the gene-dense region containing REL, the ImmunoChip study reported a long non-coding RNA, LINC01185, overlapping the lead-SNP, rs842639. However, since greater than 80% of the genome shows evidence of transcription [65], the probability of incidental overlap with long non-coding RNA is high. REL, however, is an essential transcription factor for lymphocyte development [66] and plays a critical role in autoimmune inflammation [67]. Hence, gene prioritization through integrative analyses offers not only to streamline loci discovery but also subsequent causal gene identification.
In this work, we developed a framework to predict the probability that each protein-coding gene is associated with each of 29 complex diseases. Our predictions draw on a diverse set of pathogenically-relevant relationships encoded in a heterogeneous network. The predictions successfully prioritized associations hidden from the network. Using MS as a representative example, we were able to combine our predictions with statistical evidence of association to increase study power and identify three novel susceptibility genes in this disease. The disease-specific performance (measured by the AUROC) for MS was exceeded by twelve other diseases suggesting that our predictions have broad applicability for prioritizing genetic association analyses. Prioritization can range from a genome-wide scale to a single loci where this approach can highlight the causal gene from several candidates within the same association block. For researchers focused on a specific disease, these predictions can be used to propose genes for experimental investigation. Inversely, researchers focused on a specific gene can use this resource to find suggestions for relevant complex disease phenotypes. Most previous explorations of the factors underlying pathogenicity have focused on a single domain such as tissue-specificity [68], protein interactions [69], pathways [7], or disease similarity [70]. The method presented here integrates disparate data sources, learns their importance, and unifies them under a common framework enabling comparison. Therefore, we can conclude that perturbation gene sets—the core of our top-performing feature—are an underutilized resource for disease-associated gene prioritization. Not only did perturbations encompass other set-based gene categorizations, but they greatly outperformed features based on protein interactions, pathways, and tissue-specificity, which form the basis of several prominent prioritization techniques. In addition to characterizing the overall importance of each feature, our online prediction browser visually decomposes an individual prediction into its components. We observed a prominent influence of pleiotropy, consistent with previous studies that identified pervasive overlap of susceptibility loci across complex diseases [71], especially those of autoimmune nature [72]. Since many existing prioritization techniques are agnostic to the compendia of GWAS associations, they fail to adequately leverage pleiotropy. Unlike approaches initiated from a user-provided gene list, our study only provides predictions for 29 diseases. By not relying on user-provided input, our predictions can serve as independent priors for future analyses. By predicting probabilities, we provide an extensible and interpretable assessment of association that circumvents the limitations inherent to frequentist analyses [73]. Many approaches return no assessment for the majority of genes which fall outside of their set of predicted positives. Here, we overcome this issue and provide a comprehensive and genome-wide output by returning a probability of association for each protein-coding gene. High-throughput biological data is frequently noisy and incomplete [74]. Combining orthogonal resources can help overcome these issues. Accordingly, we found that our integrative model outperformed any individual domain. While this method has shown encouraging performance, some limitations are worth noticing. For example, many biological networks preferentially cover well-studied vicinities [75]. Knowledge biases that span multiple, presumably-orthogonal resources could diminish the benefits of integration. Here, several of the literature-derived domains were removed by the lasso, suggesting redundancy. In addition, biases in network completeness can lead to high-quality predictions for well-studied vicinities and low-quality predictions for poorly-studied vicinities. The permutation analysis provided evidence of this disparity: edge-specificity was critical for top predictions yet only moderately beneficial for the remainder. Subsequently, we caution users to avoid overinterpreting predictions for poorly-characterized genes. To help place predictions in context, the online browser provides a gene’s mean prediction across all diseases and a disease’s mean prediction across all genes. However, we recognize that false negatives will continue to persist in our predictions, and users should be mindful of this limitation when interpreting results. As more systematic and unbiased resources become available [74], high-quality predictions will emerge for more network vicinities. We reason that the desirable qualities of our predictions are the consequence of the heterogenous network edge prediction methodology. The approach is versatile (most biological phenomena are decomposable into entities connected by relationships), scalable (no theoretical limit to metagraph complexity or graph size), and efficient (low marginal cost to including an additional network component). We have extended the previous metapath-based framework set forth by PathPredict [34], by: 1) incorporating regularization allowing coefficient estimation for more features without overfitting; 2) designing a framework for predicting a metaedge that is included in the network; 3) developing an improved metric for assessing path specificity; and 4) implementing a degree-preserving permutation. Metapath-based heterogeneous network edge prediction provides a powerful new platform for bioinformatic discovery.
This study was approved by the UCSF institutional review board on human subjects under protocol #10–00104. We created a general framework and open source software package for representing heterogeneous networks. Like traditional graphs, heterogeneous networks consist of nodes connected by edges, except that an additional meta layer defines type. Node type signifies the kind of entity encoded, whereas edge type signifies the kind of relationship encoded. Edge types are comprised of a source node type, target node type, kind (to differentiate between multiple edge types connecting the same node types), and direction (allowing for both directed and undirected edge types). The user defines these types and annotates each node and edge, upon creation, with its corresponding type. The meta layer itself can be represented as a graph consisting of node types connected by edge types. When referring to this graph of types, we use the prefix ‘meta’. Metagraphs—called schemas in previous work [34,35]—consist of metanodes connected by metaedges. In a heterogeneous network, each path, a series of edges with common intermediary nodes, corresponds to a metapath representing the type of path. A path’s metapath is the series of metaedges corresponding to that path’s edges. The possible metapaths within a heterogeneous network can be enumerated by traversing the metagraph. We implemented this framework as an object-oriented data structure in python and named the resulting package hetio. Users are free to browse, use, or contribute to the software, through the online repository (https: //github. com/dhimmel/hetio). The included resources, and hence the metaedges and metanodes composing our network, were selected empirically based on a balance among the following properties: 1) quality—relevance to human pathogenesis; high accuracy and an optimal trade-off between false positives and false negatives. In some cases, quality concerns prevented the inclusion of a desired metaedge. For example, we omitted ontology-based disease similarly [76] due to an inaccurate Disease Ontology hierarchy [43], and we omitted disease comorbidity due to high measurement error for uncommon diseases [77]. For included metaedges, we attempted to select the highest quality resource in that domain. 2) reusability—easily retrievable and parsable; mapped to controlled vocabularies; well documented; amenable to reproducible (scripted) analysis; free of prohibitive reuse stipulations. 3) throughput—broad domain-specific coverage generated using systematic platforms that minimize bias. While genetic interactions have previously proven informative [31], their sparse characterization in humans was deemed unfavorable for our approach. 4) diversified, multiscale portrayal of biology—capturing, in aggregate, many aspects of pathophysiology across multiple levels of biological complexity. Levels of the hierarchical architecture of biological complexity include the genome, transcriptome, proteome, interactome, metabolome, cell and tissue organization, and phenome. Balancing these considerations, we integrated as many resources as possible within our computational runtime constraints. PathPredict relied on basic logistic regression to predict coauthorship status from features corresponding to nine distinct metapaths [34]. However, faced with fewer positives to train our model and a large number of features, we adopted a regularized approach, which aims to contain the overfitting tendencies inherent to regression. Regularization penalizes complexity, a trademark of overfitting. We chose the elastic net technique of regularization [54], which is efficiently implemented for logistic regression by the R glmnet package [80]. Regularized logistic regression requires a parameter, λ, setting the strength of regularization. We optimized λ separately for each model fit. Using 10-fold cross-validation and the “one-standard-error” rule to choose the optimal λ from deviance, we adopted a conservative approach designed to prevent overfitting [80]. On the training set of gene-disease pairs, we optimized the elastic net mixing parameter (α), the DWPC damping exponent (w), and two edge inclusion thresholds. First, we optimized α and w on the 20 features whose metapaths did not include threshold-dependent metaedges. For each combination of α and w, we calculated average testing AUROC using 20-fold cross-validation repeated for 10 randomized partitionings. After setting α and w, we jointly optimized the two edge-inclusion thresholds using the AUROC for the GeTlD feature, whose metapath is composed from the two edges requiring thresholds (S3 Fig). We adopt standardized coefficients as a measure of feature effect size. Standardized coefficients refer to the coefficients from logistic regression when all features have been transformed to z-scores. Standardization provides a common scale to assess feature effect, both within and across models [81]. Starting from the complete network, a permuted network was created by swapping edges separately for each metaedge. Edge swaps were performed by switching the target nodes for two randomly selecting edges [82]. For each metaedge, the number of attempted swaps was ten times the corresponding edge count. We adopted a Markov Chain strategy where additional rounds of permutation were initiated from the most-recently permuted network [82]. A training network was generated from the first permuted network by masking 25% of the associations for testing. Testing performance for the permuted training network model is shown in S4 Fig. When contrasting this performance with the unpermuted-network model, we employed the Condensed-ROC curve to magnify the importance of top predictions [83]. Using the exponential transformation with a magnification factor of 460—the value which maps a FPR of 0. 01 to 0. 99—we concentrated on the top 1% of predictions (S4 Fig). A one-sided unpaired DeLong test [84] was used to assess whether feature-specific testing AUROCs from the unpermuted network exceeded those from the first permuted network (S2 Table). We performed a subsampling analysis for 15 gene sets—the 14 MSigDB gene sets and tissues—to assess the effect of sparsity on feature-specific performance (S7 Fig). Two without-replacement subsampling schemes were investigated: node masking and edge masking. For a specific gene set and scheme, we masked a percentage of the gene set and calculated the corresponding feature’s AUROC. We evaluated a range of percentages and performed ten subsampling repetitions for each percentage. We excluded 588 genes from the discovery phase of the multiple sclerosis analysis. First we excluded genes in the extended MHC region (spanning from SCGN to SYNGAP1 on chromosome 6 [85]) due to the complex pattern of linkage characterizing this region containing several highly-penetrant MS-risk alleles [57]. Second, we excluded putative MS genes: high-confidence primary genes from the GWAS Catalog and reported genes for the WTCCC2-replicated loci. We omitted genes in linkage disequilibrium with the putative genes by excluding: 1) consecutive sequences of nominally significant genes (using the WTCCC2-VEGAS p-values) that included a putative gene; and 2) high-confidence secondary genes from the GWAS catalog. Post exclusion, 1211 genes were nominally significant in Meta2. 5, four of which exceeded the network-based discovery threshold. Using a hypergeometric test for overrepresentation, we calculated the probability of randomly selecting 4 of the 1211 genes and Bonferroni validating at least 3 of the 4 on WTCCC2 (S3 Table). See S1–S12 Datasets for the supporting data and S13 Data for vector figures. The website provides additional resources (http: //het. io/disease-genes/downloads/) as well as an interface for browsing results (http: //het. io/disease-genes/browse/). Project related code is available from the github repository (https: //github. com/dhimmel/hetio).
For complex human diseases, identifying the genes harboring susceptibility variants has taken on medical importance. Disease-associated genes provide clues for elucidating disease etiology, predicting disease risk, and highlighting therapeutic targets. Here, we develop a method to predict whether a given gene and disease are associated. To capture the multitude of biological entities underlying pathogenesis, we constructed a heterogeneous network, containing multiple node and edge types. We built on a technique developed for social network analysis, which embraces disparate sources of data to make predictions from heterogeneous networks. Using the compendium of associations from genome-wide studies, we learned the influential mechanisms underlying pathogenesis. Our findings provide a novel perspective about the existence of pervasive pleiotropy across complex diseases. Furthermore, we suggest transcriptional signatures of perturbations are an underutilized resource amongst prioritization approaches. For multiple sclerosis, we demonstrated our ability to prioritize future studies and discover novel susceptibility genes. Researchers can use these predictions to increase the statistical power of their studies, to suggest the causal genes from a set of candidates, or to generate evidence-based experimental hypothesis.
Heterogeneous Network Edge Prediction: A Data Integration Approach to Prioritize Disease-Associated Genes