key: cord-0023089-iaqzms3g authors: Buphamalai, Pisanu; Kokotovic, Tomislav; Nagy, Vanja; Menche, Jörg title: Network analysis reveals rare disease signatures across multiple levels of biological organization date: 2021-11-09 journal: Nat Commun DOI: 10.1038/s41467-021-26674-1 sha: 58b7296eb9f02cd30a4ab97e9eec6d6bf376cba9 doc_id: 23089 cord_uid: iaqzms3g Rare genetic diseases are typically caused by a single gene defect. Despite this clear causal relationship between genotype and phenotype, identifying the pathobiological mechanisms at various levels of biological organization remains a practical and conceptual challenge. Here, we introduce a network approach for evaluating the impact of rare gene defects across biological scales. We construct a multiplex network consisting of over 20 million gene relationships that are organized into 46 network layers spanning six major biological scales between genotype and phenotype. A comprehensive analysis of 3,771 rare diseases reveals distinct phenotypic modules within individual layers. These modules can be exploited to mechanistically dissect the impact of gene defects and accurately predict rare disease gene candidates. Our results show that the disease module formalism can be applied to rare diseases and generalized beyond physical interaction networks. These findings open up new venues to apply network-based tools for cross-scale data integration. O ver the past 2 decades, rapid advances in DNA sequencing technology allowed us to uncover the genetic basis of over 6000 rare diseases [1] [2] [3] . In contrast to common diseases, which are typically characterized by a complex interplay between multiple genetic and environmental factors, rare diseases can often be pinpointed to a single genetic lesion. Rare diseases thus offer unique opportunities to mechanistically dissect the relationship between genetic aberrations and their phenotypic consequences, which can then inform targeted treatment strategies. For individual rare diseases, this potential for a molecularly rooted, personalized medicine could already be demonstrated, for example in rare immunodeficiencies 4-6 , neurodevelopmental 7, 8 , and metabolic disorders 9, 10 . At the same time, the costs and extended timelines of these individual efforts also highlight the need for novel, systematic approaches for investigating the large number of rare diseases that still remain uncharacterized. To this end, several practical and conceptual challenges need to be overcome: First, rare disease phenomena cover a wide spectrum, from highly cell-type or organ-specific phenotypes to heterogeneous, syndromic diseases that affect the whole body. Our understanding of how a genetic aberration impacts various scales of biological organization between genotype and clinical phenotype is very limited. Second, the enormous complexity within and between different organizational scales, such as the transcriptome, proteome, intra-or intercellular communication, also poses important technical challenges: How can we identify and integrate the most relevant data? Third, the rarity of many conditions with monogenic origins implies that data are usually scarce. Traditionally, rare diseases have been studied following a one-gene, one-pathway, one-disease paradigm. A systematic approach for transferring knowledge from one rare disease to another, and for investigating differences and commonalities between different diseases, is still missing. In this work, we propose a network-based framework for systematically investigating rare diseases that addresses these challenges, and, in turn, use the large number of rare diseases with a well-described genetic origin to deepen our understanding of disease-associated perturbations of molecular networks. Specifically, we introduce a multiplex network approach for integrating different network layers that represent different scales of biological organization ranging from the genome to the transcriptome and the phenome. A systematic characterization of the network signatures of all rare diseases with known genetic causes allowed us to identify the connectivity patterns that determine the importance of a particular scale of biological organization for a given rare disease. Finally, we explored how these systems-level insights may help contextualize individual genetic lesions, investigate the impact of disease heterogeneity, and be translated into clinically actionable tools for the genetic diagnosis of rare disease patients with unknown gene defects. Constructing a gene network bridging molecular and phenotypic scales. Rare diseases affect many scales of biological organization which, conversely, may provide valuable information for elucidating a particular gene defect. At the genetic level, for example, interplay between genetic variants can modulate phenotypic outcomes 11 or even completely rescue disease-associated variants 12 . At the protein level, members of the same complex or pathway are often implicated in similar phenotypes 13, 14 and expression patterns of a particular gene can reveal affected cell types and tissues [15] [16] [17] . Finally, phenotypic similarities with known human or animal model gene defects can guide the annotation of genetic variants with unknown consequences 18 . To integrate these diverse relationships into a unifying, genecentric framework, we constructed a multiplex network comprised of several layers: The nodes in each layer represent genes, the links represent their respective relationship at a particular scale of biological organization, ranging from direct interactions between gene products at the molecular level to phenotypic similarity of associated diseases at the phenotype level ( Fig. 1a, b) . We compiled information from seven databases and, where appropriate, applied a range of techniques for extracting gene relationship, such as bipartite mapping, ontology-based semantic similarity metrics and correlation-based relationship quantification, as well as filtering based on both statistical and network structural criteria 19 (Fig. 1c, d 24, 25 . Characterizing the network architectures across biological scales. To characterize the resulting cross-scale gene relationships, we first quantified the global similarity between all pairs of network layers A and B by the overlap of their respective sets of edges E: S AB ¼ jE A \ E B j=minðjE A j; jE B jÞ. The highest similarities were found within the transcriptomic scale: co-expression networks of different tissues have an overlap of up to S ¼ 0:49 (between brain tissues), compared to an average similarity of S ¼ 0:05 between networks of other scales. A major contribution to this elevated similarity is given by a core of links that is preserved across multiple tissues. We found that the proportion of links that connect essential genes increases with the number of tissues in which a particular link is present (Fig. 1e ). This suggests that the common core is related to essential housekeeping activities. To represent pan-tissue and tissue-specific interactions separately, we extracted broadly preserved co-expression edges and considered them as a separate core transcription network layer, consisting of 12,364 nodes and 1,062,924 edges ( Supplementary Fig. 1c , d, Methods). We further combined redundant tissue types, resulting in a final set of 38 tissue-specific networks used in the downstream analyses (Fig. 1f, Methods, Supplementary Data 3) . These tissue-specific networks still form a recognizable cluster within the multidimensional scaling (MDS) projection of the relative similarities between all networks (Fig. 1g ). The differences between tissues, however, are comparable with differences to networks of other scales (median similarity among tissues: S ¼ 0:043; similarity to other scales: S ¼ 0:018). The clear separation between most network layers (median similarity S ¼ 0:033) indicates that each layer contains unique information ( Supplementary Fig. 3b ). At the same time, a comparison with randomized networks reveals that a significant amount of interactions are preserved across levels of organization (Supplementary Fig. 3c, d) , as shown by a significant similarity for 96.5% of all network pairs (empirical p-value < 0.05, see Methods). Finally, we noticed that the relative position of all network layers in Fig. 1g suggests a representative role for the protein-protein interaction (PPI) layer, which is located in a central position and close to the layers that directly encode phenotypic similarities. We next compared the networks at the different biological scales in terms of five structural characteristics: genome coverage, overall connectivity, clustering, assortativity and literature bias (Fig. 1h, Methods) . The results revealed a wide structural diversity: The network layer with the highest genome coverage is the PPI scale, covering 17,944 proteins. This is due to the combination of a large number of literature curated small-scale experiments and several large-scale screening efforts. Such systematic, genome-wide measurements also underlie the high coverage of the transcriptomic layers (with a total number of N = 17,432 genes across all tissues, and an average number of 10,527 genes per tissue, Supplementary Fig. 1d , see Methods for the filtering processes). Our incomplete understanding of how these molecular interactions translate into biological processes, however, is indicated by the low coverages observed among the functional and phenotypic levels (N = 2407 and 3342 for the molecular function and HPO networks, respectively). The high connectivity and clustering among these functional layers, in turn, is the basis for their predictive power for transferring gene annotations within functional clusters 11,20 (e.g., edge density = 1:13 10 À2 and clustering = 0.73 for the co-essentiality network). The PPI represents the sparsest network (edge density = 2:359 10 À3 ; average density across all layers = 7:76 10 À3 ), which, in part, reflects the incompleteness of currently available data 26 . Curiously, the PPI is the only network in our collection that exhibits a (modest) level of disassortativity (r ¼ À0:08), i.e., a tendency of hubs to connect preferentially to low-degree nodes, a property that was previously suggested to be a universal feature of biological networks 27 . Disassortativity may arise when the neighbors of high-interest nodes are mapped out more extensively than the interaction partners of these neighbors. For the PPI, this is likely to be the case in network data curated from hypothesis-driven, small-scale experiments, but can also occur in unbiased large-scale efforts ( Supplementary Fig. 3e , f, Methods). A further characterization of curated and unbiased subsets of the PPI ( Supplementary Fig. 4a , b, Methods) revealed that the relatively high literature bias present in the PPI, as measured by the correlation between the degree of a protein and the number of associated publications (Spearman's ρ = 0.59, Supplementary Fig. 4b ), is largely driven by its curated subset, which represents 87% of the full PPI. This emphasizes that despite recent efforts towards unbiased, high-throughput protein-protein interaction screening, a large fraction of the currently available PPI network information still reflects the reductionist, hypothesis driven research paradigm where new knowledge preferentially accumulates around proteins with an already known important function. This literature bias is notably absent in all other network layers. In summary, the structural diversity observed among the individual network layers reflects both organizational principles intrinsic to a particular biological scale, as well as technical or historical details pertaining to the curation process of the underlying database ( Supplementary Fig. 4c-f ). We expect that this diversity further corresponds to complementary pieces of information contained in the different biological scales, collectively increasing their potential to drive novel insights into the relationships between rare disease genes. Identifying cross-scale network signatures of rare diseases. To investigate the connectivity patterns among rare disease genes, we collected 3953 genes associated with 3771 rare disease terms from the Orphanet database, the largest rare disease ontology and resource for genetic associations (Supplementary Data 4). Collectively, rare diseases represent an extraordinarily rich resource of causative genetic aberrations and their phenotypic consequences. For individual rare diseases, however, the situation is the opposite: Over 3501 diseases in the Orphanet database (~93%) are associated with fewer than five genes. This represents a major challenge for systematic, comparative rare disease research in general, and for network-based approaches in particular: Network approaches are based on the fundamental observation that genes associated with the same disease are not scattered randomly in molecular networks, but aggregate in disease-specific neighborhoods or "disease modules" 26, 28 . However, the incompleteness of currently available network maps sets a lower bound for the number of genes that can be recognized as a connected module. This minimal number was estimated to be around 20 for the PPI network 26 , so that individually, only few rare diseases have a sufficiently large number of associated genes. We hypothesized that the disease module concept can be generalized to groups of rare diseases with closely related phenotypes. Collectively, these related rare diseases could thus reach the required minimum number of genes to form a recognizable disease module (Fig. 2a) . To test this hypothesis, we used the hierarchical classification of rare diseases within the Orphanet Disease Ontology to aggregate rare diseases with similar phenotypes and collect all genes associated with their corresponding descendant terms. We identified a total of 26 rare genetic disease groups that are sufficiently broad or well-studied, respectively, to result in a number of associated genes required for network module approaches (i.e., more than 20), while retaining the pathophysiological specificity of rare disease phenotypes ( Supplementary Fig. 5a ). The disease groups range from smaller groups, such as RASopathy (ORPHA:536391) or rare genetic vascular diseases (ORPHA:233655) (with 20 and 22 associated genes, respectively), to large groups with over 1000 associated genes, such as rare genetic neurological disorder (ORPHA:71859) or rare developmental defect during embryogenesis (with 1649 and 1598 associated genes each). The average number of genes per disease group was 339 ( Fig. 2b and Supplementary Data 5). Despite the wide range in the total number of associated genes per disease group, the average number of genes per disease term remains comparable across all disease groups, thus ensuring similar levels of disease specificity across the disease domain. In addition, there is only little overlap between the disease terms contained in the different groups, with 90.5% of all disease pairs being distinct (Jaccard Index < 0.1), indicating that the groups provide non-redundant disease definitions ( Supplementary Fig. 5b) . We first inspected the network localization of the aggregated rare disease groups within two-dimensional network embeddings obtained from the node2vec algorithm 29 , which aims to preserve The multiplex network consists of 46 network layers, each representing a particular type of gene relationships, ranging from genetic interactions to phenotypic similarity. c Methods used for inferring networks: bipartite mapping was used to build gene relationships based on common annotations, e.g., pathways; semantic similarity was used to define relationships based on annotation similarity; correlation analyses were used to identify co-expression. d Weighted and dense networks were subsequently filtered based on structural network criteria for extracting the most relevant interactions. e Co-expressed gene pairs found in a higher number of tissues tend to be essential, reflecting core cellular functions. Edges found in five or fewer tissues were considered tissue-specific. f Full co-expression profiles are highly similar between tissues and thus redundant (lower triangle). The removal of core transcription profiles reveals tissue-specific patterns of the co-expression networks (upper triangle). g The multi-dimensional scaling (MDS) plot based on edge overlap similarity of all networks shows a clear distinction between major types and subtypes of the network layers. h Major network characteristics for all considered network layers: number of nodes edge density, global clustering, assortativity and social bias, as measured by the correlation between node degree and number of associated publications. The values of the 38 individual co-expression networks are shown in the form of a distribution. network distances between nodes ( Supplementary Fig. 5c , Methods). Fig. 2c shows the resulting network landscape of the human phenotypic network with four rare disease groups highlighted: rare bone, immune, hematologic and neurological diseases (the complete landscape of all diseases in all networks can be explored via our MultiOmeExplorer web app: www.menchelab.com/MultiOmeExplorer, see also Supplementary Fig. 6 ). We found that all four disease groups localize within specific network neighborhoods. Given that the HPO network is based on phenotypic similarity of individual gene defects, this localization confirms that aggregating diseases based on disease ontology relationships indeed leads to groups of phenotypically related diseases. We further noticed that the different disease groups cover network areas of varying size, from highly localized immune diseases to more broadly spread neurological diseases. In part, this spread can be attributed to the larger number of genes associated with the latter disease group. More generally, it may reflect varying degrees of coherence and specificity among the phenotypic manifestations of the diseases represented within a particular group. The close relationship between the spread of disease-associated perturbations within molecular networks and the heterogeneity of clinical symptoms has previously been shown for complex diseases 30 . Similarly, the spread of the rare neurological disease cluster recapitulates the high level of comorbidities observed among affected patients. Finally, we noted that the proximity between neighborhoods is indicative of disease similarity, e.g., between rare immune and hematologic diseases, where the interplay between blood and immune system often leads to similar phenotypes. We next inspected the network signatures of rare disease groups across different network layers. Figure 2d shows that rare genetic cardiac diseases are strongly localized on a heart-specific co-expression network (heart right ventricle; HRV) and the human phenotypic similarity network (HPO). The more dispersed signals on the PPI network and the network of shared biological processes (GO:BP), on the other hand, suggest that the respective genes might be involved in a broad range of molecular processes that cannot be adequately depicted in a twodimensional projection. Quantifying network modularity of rare diseases. The results so far indicate that the concept of disease modules, observed widely across complex diseases on PPI networks, can also be generalized to groups of rare diseases and to other network data representing relationships beyond the molecular scale of PPIs. Based on the heterogeneous degrees of modularity for different diseases and networks observed above, we further hypothesized that the degree of modularity can be related to the degree of relevance of the underlying information to a particular disease phenotype. To investigate this hypothesis and further dissect the characteristics of rare disease modules across biological scales, we systematically assessed all rare disease groups across all network layers. We quantified the level of modularity by the significance of the size of the largest connected component (LCC) of disease genes on a given network, as measured by the corresponding z-score compared to random gene sets (Fig. 3a , Methods). Figure 3b shows the module significance for all rare disease groups on all network layers summarized in one heatmap. We observed a high degree of differential modularity, i.e., the levels of localization vary greatly between disease groups and network layers. The largest number of significantly localized rare disease groups are found on the PPI network, the phenotypic networks (HP, MP), the core transcription network, and the network of shared biological processes. This consistent localization across a wide range of rare diseases confirms the existence of disease modules also for rare diseases. In contrast to the core transcription layer observed to be relevant across multiple disease groups, the tissue-specific co-expression networks provide a more disease-specific picture with unique signatures that reflect the molecular mechanisms that underlie a particular disease group on a given tissue. For example, the wider localization pattern of rare neurological disease genes in the Rare disease grouping and network mapping reveal network-and disease-specific connectivity patterns. a Rare genetic diseases are typically associated with only a few genes and therefore remain fragmented on molecular networks. Grouping rare diseases by phenotypic similarity can overcome data scarcity and result in identifiable disease modules, thus allowing for further network-based inspection. b Voronoi treemap showing the 26 rare genetic disease groups used in this study. The size of each disease group is proportional to the number of associated genes. c Network landscape obtained using the node2vec embedding algorithm. Network distances between genes are preserved in the embedding and illustrate differential modularity of different rare disease groups on the Human phenotype similarity network layer. The bright dots represent disease associated genes and the blue contour map represents all genes in a network. d Localization of the rare cardiac disease group on different network layers. Fig. 2c corresponds to their significant modularity across co-expression networks in a wide range of tissues, which, in turn, reflects the often syndromic and heterogeneous phenotypes of these diseases. Using differential modularity to contextualize rare disease gene clusters. The individual layers within the cross-scale network capture different pathobiological mechanisms. The observed differential modularities can thus offer insights into the disease etiology specific to a particular layer. For example, rare genetic gastroenterological diseases, a disease group consisting of 92 disease terms with 140 associated genes, were found to be significantly localized on five network layers (Fig. 4) . Detailed inspection and enrichment of these submodules ( Supplementary Fig. 6a , b, Methods) enables us to interpret the disease characteristics within each layer: We found that genes causing the Bardet-Biedl syndrome (BBS) form pronounced clusters in the phenotypic and PPI layers. Together with the absence of modularity in other layers, this pinpoints that the emergence of this particular disease phenotype is mainly determined by interactions at the protein level, while co-essential, functional or pathway levels play less important roles. This observation is supported by our current knowledge of BBS pathological mechanisms: The proteins encoded by BBS genes form a complex crucial for transporting vesicles to cilia, a process whose defect is suspected to be a major cause of BBS 31 . At the same time, these proteins are of diverse functional character 32 and involved in disparate pathways 33 , explaining the lower modularity on the respective network layers. shows the modularity of all rare disease groups across all network layers as measured by the respective module size significance (p < 0.05: *, p < 0.01: **, p < 0.001: ***; p < 0.0001: ****, Benjamini-Hochberg corrected empirical p-values determined by node randomization, see Methods). In the tissue-specific network layers, only selected disease groups display pronounced modularity, often recapitulating known mechanisms and tissue specificities of particular rare diseases, but also revealing novel relationships. Network layers containing relationships that are relevant across biological levels of organization, such as protein-protein interaction, phenotypic and functional similarity networks, also display modularity across a wide range of disease groups. We also observed differential modularity among rare cancerrelated disorders. Pancreatic carcinoma genes are involved in various apoptotic processes and significantly connected on both protein and pathway levels, highlighting the BAX and survivin protein complexes, and FGFR1, SCF-KIT, and PGDF signaling pathways, respectively. Modulation of gene expression involved in these processes has been reported to play a key role in pancreatic tumor growth 34, 35 . Moreover, genes involved in cancers of gastrointestinal tracts (esophageal, gastric and colorectal cancer) form a strong cluster on the esophagus co-expression network, revealing their interconnected roles at the gene regulatory level. Another notable cluster among rare gastroenterological diseases is given by the Shwachman-Diamond syndrome: 90% of all patients have a mutation in the SBDS gene, which is known to be involved in ribosome maturation, but otherwise poorly annotated. Recent studies indicate that the remaining 10% of patients with similar disease phenotypes have mutations in other genes involved in ribosome biogenesis 36, 37 . The interactions between these genes and SBDS are only captured in the genetic (co-essentiality) and the transcriptomic (co-expression) layers. These observations pave the way for a detailed mechanistic interpretation of how different network layers contribute to the etiology of individual rare diseases, as well as for identifying mechanisms that are shared among phenotypically related rare diseases. Modularity as quantification of pathobiological relevance. Our findings suggest that the significance of the network localization of rare disease genes on a particular layer of the cross-scale network may be used to quantify the pathobiological relevance of the respective level of biological organization for the disease. The information in Fig. 3b could thus be interpreted as a networkdisease relevance score (π) for distinguishing information that truly reflects the system of interest from unrelated information and potential noise. Specifically, we hypothesized that if the associated genes of a disease group are significantly connected on a particular network, then the network has predictive power for discovering novel genes associated with the disease. In contrast, if the genes are scattered on a particular network, then it is likely uninformative for the discovery process. To test his hypothesis, we developed an informed multiplex network propagation algorithm, in which the overall probability p m to visit a given layer m out of all L layers, p m ¼ ∑ L i¼1 pðijmÞ, is proportional to its respective level of relevance π m , which can be achieved by incorporating the detailed balance condition π m pðmjnÞ ¼ π n pðnjmÞ (Methods). To validate the potential of this informed propagation algorithm for disease gene discovery, we performed a 10-fold cross-validation for the retrieval of associated genes for all disease groups and assessed the performance through the area under the receiver operating characteristic curve (AUROC, Fig. 5a) . We compared four different scenarios incorporating (i) only the single most informative network (i.e., the one with the highest number of significantly localized disease groups; here: the PPI), (ii) the single most relevant network for each disease group ( Supplementary Fig. 7a) , (iii) all networks and (iv) only the most relevant networks (i.e., networks with a modularity significance of p-value < 0.05, Benjamini-Hochberg correction for multiple hypotheses). We found that all four different sets of networks performed reasonably well, with AUROC ranging from 0.65 to Fig. 3) , which capture relevant relationships on different scales. Diseases in the disease group exhibit unique connectivity patterns in the network layers where their disease characteristics can be derived. For example, a strong phenotypic cluster of Bardet-Biedl syndrome genes can be derived from the protein complex (BBSome) whose defects lead to cilia dysfunction, while pancreatic carcinoma cluster is most observable on the pathway level where its causal genes interact physically and are also member of crucial signaling pathways. The treemap represents disease entries in Orphanet (leaf terms with at least two gene associations) that belong to the rare genetic gastroenterological disease group (69 diseases with only one gene association are not shown separately). (Fig. 5b) , confirming the general applicability of multiplex network propagation to rare disease gene prediction. A comparison between the four methods revealed that incorporating only relevant network layers (median AUROC = 0.90) generally outperforms the PPI (median AUROC = 0.73), the most relevant single layer benchmark (median AUROC = 0.79), as well as the incorporation of all layers (median AUROC = 0.86), with corresponding Bonferroni-Holm corrected Durbin-Conover test p-value = 3e-16, 1.22e-6, and 0.002 respectively (Fig. 5c) . We concluded that network modularity thus provides a network-based criterion to curate and integrate the most relevant data and levels of biological organization for a specific disease. Interestingly, we also observed differences in retrieval performance related to characteristics of the diseases themselves: We found that syndromic disease groups, i.e., those with significant disease modules across multiple tissue types, tend to have lower retrieval performance and benefit from incorporating all tissue co-expression networks (Spearman's ρ = −0.53, p-value = 0.004, Fig. 5d left panel, Supplementary Fig. 7b ). On the one hand, this poses a challenge for disease groups that manifest various anatomical features such as rare genetic neurological disorders. On the other hand, this reflects limitations of broad ORDO disease group definitions such as rare genetic defects during embryological development. We further found that the retrieval performance correlated negatively with the number of genes associated with a particular disease group (Spearman's ρ = −0.83, p-value = 1.94e-6, Fig. 5d right panel) . These two factors are closely related, as both the syndromicity level and overall heterogeneity tend to increase as more genes are involved in the disease group. Taken together, these findings indicate that well defined disease groups with low to moderate number of associated disease genes are more likely to capture molecular Using network modularity as relevance prior in the informed propagation algorithm for gene prioritization. a Schematic overview of the informed multiplex network propagation algorithm that incorporates modularity as measure of relevance of a particular network level for a given disease group. b Comparison of 10-fold cross-validation performance in rare disease gene retrieval for different choices of included networks: Informed algorithm with most relevant network (blue), all networks (green), the PPI (red), and the single most relevant layer for each disease (yellow). Dashed lines show median value across all folds, shaded areas represent the interquartile range. The retrieval performance indicates that disease mechanisms are generally better recapitulated by incorporating relevant networks only. c Comparison of the AUROC from all four methods. Utilizing the significant networks lead to more accurate disease gene retrieval compared to all networks, the single most relevant layer, or the PPI. (Bonferroni-Holm corrected Durbin-Conover test pvalue = 0.026, 1.22e−6, and 3e-16 respectively). Threshold for p-values: p < 0.05:*, p < 0.01:**, p < 0.001:***, p < 0.0001:****; n = 26 rare disease groups across all network sets. Bounds of box represent 25th and 75th percentiles, center the median, whiskers 10th and 90th percentiles, respectively. d Factors correlated with the retrieval performance. The algorithm that incorporates all networks can outperform the informed algorithm for diseases with high levels of syndromicity, i.e., disease that manifest in multiple physiological systems (left, Spearman's ρ = −0.53, corresponding p-value = 0.004). Decreasing functional relevance as the number of genes increases also led to lower predictive performance (right, Spearman's ρ = −0.83, p-value = 1.94e-6). The corresponding p-value of correlation was determined by Fisher z-transformation, two-sided. disease characteristics at a level of specificity that results in better network-based predictions. This suggests that more fine-grained, mechanism-based disease definitions, together with highresolution phenotyping will aid in further improving the predictive power of the introduced network methods. To further dissect the contribution of individual networks and potential curation biases on the overall predictive power, we performed several additional benchmarks on different subsets of the multiplex network (Methods). Our comparisons between curated, unbiased and size-matched random subsets of the PPI indicate that the performance is largely driven by network size rather than potential literature biases in the interaction curation process (Supplementary Fig. 7c) . We also evaluated the differences in performance upon removing individual layers, as well as groups of layers from the full multiplex network (Supplementary Fig. 7d) . The results suggest that the performance is not driven by individual network layers and that the predictive power of the multiplex network can be best understood as a collective characteristic of all disease relevant layers. Application to candidate gene prioritization in rare disease patients. Based on the performance of the informed multiplex propagation for retrieving genes across all rare disease groups we hypothesized that the method can also act as an additional evaluation metric for prioritizing genomic variants in individual rare disease patients. Starting point in a diagnostic setting is nextgeneration sequencing of a patient's genome, typically yielding rare genomic variants (allele frequency < 1%) in dozens to hundreds of different coding regions, and with unknown consequences. These variants may be further filtered down, for example based on frequency in the general population, deleteriousness scoring, or segregation analysis, resulting in up to a few dozen high confidence candidate genes 38, 39 . Identifying the one causal gene among them remains a critical challenge both in research and in clinical practice (Fig. 6a) . We tailored the informed propagation algorithm to individual patients by using seed genes associated with patient-specific phenotypes, combined with the network relevance scores from the corresponding Orphanet disease group ( Supplementary Fig. 7e) . Altogether, this enables us to perform patient-specific multiplex network propagation to prioritize candidate genes. We applied the method to filtered lists of genes with rare variants obtained from a cohort of 139 rare disease patients suffering from various neurological symptoms with intellectual disability as a predominant phenotype (Fig. 6b, Supplementary Fig. 7f , g, Supplementary Data 7, and Methods for details of the cohort). The causal variants of all patients were already confirmed and could thus be utilized for benchmarking. After standard methods of filtering for high confidence variants were exhausted, up to 997 candidate genes per patient remained (mean = 401.2). We found that our algorithm prioritizes causal genes with an overall AUROC of 0.95 (Fig. 6c) . Furthermore, we benchmarked the performance of our method against predictions based on various gene-level properties, using the same data as in the network construction, specifically (1) pathway information, (2) expression level information, (3) literature counts and (4) phenotypic similarity (see Methods for details). Among these methods, phenotypic similarity was most predictive (AUC = 0.87), followed by literature counts (AUC = 0.72), expression information (AUC = 0.71) and pathway counts (AUC = 0.59). However, the informed multiplex propagation outperformed all gene-based methods (p-value = 8.39e-5, DeLong's test of ROCs between our approach and the best performing gene-level method, i.e., phenotypic similarity). Note that for the specific use case of gene variant prioritization in rare disease patients, it is instructive to not only consider the global ROC-based performance, but also the exact ranking of the causal gene when assessing the utility of the framework in research and clinical settings. In our cohort, the multiplex propagation method placed the true causal gene among the top five ranked genes for 64 out of 131 patients (48.9%). For the purely gene-based methods, the causal gene was among the top five in only between 4 and 11 patients (3.1-8.4%, Fig. 6d) . We also performed an additional benchmark using a temporalholdout setting to ensure that the performance of our method is not primarily driven by confirmatory biases ( Supplementary Fig. 8a, Methods) . To this end, we curated a set of 21 patients with causal genes that were unknown at the time of network construction, thus minimizing the likelihood that their disease association contributed to information in any of the curated databases ( Supplementary Fig. 8a, Methods) . We found that the overall performance as measured by the AUROC remained high for all tested prediction methods. The observed slight reductions were within the 10-fold interquartile range in most cases and may also be attributed to the smaller sample size. For example, the informed multiplex propagation AUROC was reduced from 0.90 to 0.86 ( Supplementary Fig. 8b) . A closer inspection of the ranking showed that our framework maintained its proportion of true causal genes being ranked in the top five gene list, whereas almost all gene-based approaches had difficulties in retrieving them at highly ranked positions ( Supplementary Fig. 8c, d) . In the context of complex diseases, numerous network-based studies have revealed an intimate relationship between genetic disease associations, their interaction patterns, and pathophysiological manifestations 40, 41 . Most importantly, it was found that disease genes are not scattered randomly in molecular networks, but instead agglomerate in disease-specific modules 26 . Molecular networks can thus serve as maps to guide the search for new disease genes [42] [43] [44] , suggest drug repurposing [45] [46] [47] and combination strategies 48, 49 , or elucidate disease relationships 26, 50 , to name but a few important applications. Our work expands this concept in several directions: First, we showed that by aggregating individual gene defects into groups of related phenotypes, we can apply tools originally developed for common, polygenic diseases also to rare, monogenic diseases. Our comprehensive analysis of over 3,584 individual gene defects revealed that as a group, they exhibit network signatures similar to those observed for complex diseases. This opens up a wide range of network medicine tools and concepts to be applied to rare diseases. Existing tools, for example for prioritizing rare variant genes, often augmented by additional clinical data 51-55 , demonstrate the potential for network-based methods in this area. We further showed that the central network medicine concept of disease modules can be generalized towards multiplex networks representing various levels of biological organization. Previous work relied prominently on physical protein-protein interactions, which have been mapped out systematically for nearly two decades 13, 14, 56, 57 . Our analysis of 46 network layers containing over 20 million interactions showed that disease modules can be identified across a wide range of relevant gene relationships. We further found that the degree of modularity is indicative of the impact of disease-associated perturbations on a particular level of biological organization, and thereby determines the disease relevance of datasets from the respective level. The performance of the informed propagation algorithm for rare disease gene prediction demonstrates the practical utility of this finding. We expect that the general principle for identifying the most relevant datasets will be applicable in other contexts as well, including studies on cancer and other complex diseases. Indeed, as biomedical research is becoming more data intensive in general, and more biological network maps become available in particular 58, 59 , new strategies for integrating diverse data are required. A range of methodologies have been developed for this purpose, including network-based strategies 60-63 and advanced machine-learning approaches [64] [65] [66] . Our results could potentially enhance these strategies by using disease modularity as a criterion for curating and excluding potentially uninformative datasets. Finally, to enable a broad community of researchers in the areas of rare disease, network medicine or biomedical data integration to build on our work, all datasets and algorithms presented in this work are publicly available. Resources and network construction. Resources used in the multiplex network construction are listed in Table 1 . We incorporated seven major databases, each representing distinct biological layers. Protein-protein Interactions. Protein-protein interaction data was taken from the HIPPIE database 21 and filtered for interactions with supporting PubMed articles. To assess the impact of interactions collected from small-scale, hypothesis-driven experiments compared to those stemming from large-scale, unbiased screens, we further collected the most recent versions of the two largest systematic high-throughput PPI studies: the Human Reference Interactome (HuRI) based on yeast two-hybrid (Y2H) screening 13 (retrieved from http://www.interactome-atlas.org/ data/HuRI.tsv on 31 May 2021), and the BioPlex interactome constructed from affinity-purification mass spectrometry profiling 67 (retrieved from https:// bioplex.hms.harvard.edu/data/BioPlex_293T_Network_10K_Dec_2019.tsv, on 31 May 2021). The PPI network layer can therefore be split into two categories: the unbiased PPI for interactions that are contained in any of these two resources, and the curated PPI for the remaining edges ( Supplementary Fig. 4a) . Tissue Co-expression networks. Transcriptomic data is one of the most abundant publicly available high-throughput data. Differential expression profiles across tissues and cell types have been widely analyzed as a probe for disease specificity. In the context of network analyses, expression data has been used in two major ways: as a means to filter out genes from generic interactomes based on their expression level in a particular context of interest 13, 60 , and for constructing co-expression networks. Here, we follow the latter approach, and use co-expression as a proxy for tissue-or cell type-specific functional and regulatory relationships. As primary resource we used the Genotype-Tissue Expression (GTEx) data 68, 69 , which provides genome-scale expression profiles across 53 human tissues that have been used previously to construct co-expression networks 15, 16 . We used the following pipeline: 1. We downloaded the GTEx expression profiles in the format of transcripts per million (TPMs) from the Expression Atlas (https://www.ebi.ac.uk/gxa/ experiments/E-MTAB-5214/). The data was subsequently processed using the bioconductor package SummarizedExperiment (https:// bioconductor.org/packages/release/bioc/html/ SummarizedExperiment.html). 2. Tissues with a number of samples less than a minimum quality threshold similar to the GTEx Portal (n ≥ 70) were removed. These included fallopian tube (n = 14), ectocervix (n = 12), endocervix (n = 10) and urinary bladder (n = 24). 3. For the remaining 49 tissues, we further merged tissues with similar expression profiles to reduce redundancy and increase signals 70 . Most notably, brain regions (13 tissues) were merged into three major groups and relabeled by their anatomical entities (Supplementary Fig. 1a) . This process not only merged potentially redundant tissues, but also increased sample sizes for some tissue groups that would otherwise have been undersampled. Fig. 1b) . 5. For each tissue, Spearman's correlation coefficient (ρ) of all protein-coding gene pairs was used to determine the strength of their respective coexpression levels. Gene pairs with |ρ| ≤ 0.75 were discarded, resulting in 11,161 ± 1082 genes per tissue. 6. We applied a disparity filter 19 to remove weak, structurally redundant edges and to extract the backbone of each network. Edges with a corresponding disparity filter p-value < 0.05 were selected. This process yielded 10,526 ± 1825 genes per tissue. Note that even though the number of nodes decreased only slightly, the disparity filter excluded a large amount of spurious correlations (median number of interactions before and after = 1.83e6 and 4.78e5, respectively, Supplementary Fig. 1b (right) ). The disparity filter represents a dynamic cutoff, where lowly expressed genes tend to be removed and highly expressed genes tend to remain, while also allowing for the detection of lowly expressed genes that are strongly correlated with other genes (Supplementary Fig. 1e) . As a reference, we also show the comparable reduction of remaining genes if they were filtered using a standard expression threshold of TPM > 1 (13,567 ± 874 genes, Supplementary Fig. 1b) . The resulting networks consist of edges that are shared across multiple tissues (core transcriptional modules), as well as edges that are only present in a small number of tissues (tissue-specific modules). We considered edges present in less than five tissues as tissue-specific, and edges present in at least five tissues as core transcriptional modules (Supplementary Fig. 1c, d) . Ontology-derived functional and phenotypic similarity network. To capture gene relationships on functional and phenotypic levels, we incorporated expert curated data and systematic ontologies. To transform ontological annotations into genecentric networks, we defined that two genes are functionally or phenotypically connected if they are semantically similar based on the corresponding ontology 71, 72 as follows: We first compared several widely used measures of semantic similarity to ensure that the scores are robust for our purposes: 1. Information content (IC)-based similarity based on Resnik's method 73 . The similarity of two terms is derived from their most informative common ancestor (MICA) in the ontology. Given ontology terms t 1 and t 2 , their pairwise similarity is given by sim Resnik ðt 1 ; t 2 Þ ¼ ICðMICAÞ, where ICðtÞ ¼ ÀlogðpðtÞÞ, p(t) represents the frequency of term t defined by pðtÞ ¼ n t N ; n t denotes the number of descendants of term t, and N the number of descendants of the root term of interest in the ontology tree. 2. Information content (IC)-based similarity based on Lin's method 74 Resnik's method, Lin's similarity measure restricts the value to be in the range between zero and one, and is given by sim Lin ðt 1 ; t 2 Þ ¼ 2ICðMICAÞ ICðt 1 ÞþICðt 2 Þ 2 ½0; 1. 3. After collecting all pairwise term similarities for annotations of two genes, we next employed the Best-Match Average (BMA) strategy to combine them into a single gene similarity score. Their pairwise similarity of genes g 1 and g 2 with m and n annotated terms, respectively, is given by ;where S 2 R mxn is the matrix containing the pairwise similarity values of the ontology terms associated with the two genes, rowmaxðSÞ and colmaxðSÞ are vectors of length m and n, containing the maximum similarity values across all rows and columns of matrix S. 4. Frequency-based similarity, where the similarity between two genes is given by the number of shared annotations, i.e., sim freq ðg 1 ; where T gk is the set of ontology terms (including ancestor terms) associated with gene k. We found that the respective similarity values are strongly correlated, indicating that the resulting networks are robust against details of the used methods ( Supplementary Fig. 2a) . We chose to proceed with the IC-based Resnik's method with the Best-Match Average (BMA) combination strategy, as it has been demonstrated to both be among the simplest methods, while also providing the most reliable performances across different tasks 71, 75 . We used the R packages GoSemSim 76 and OntologyX 77 to navigate and compute the similarity measurements. Gene pairs with minimal similarity value, i.e., pairs whose only common annotation is the root term of the considered ontology branch (i.e., "Molecular Function" or "Biological Process") were considered as unrelated and therefore removed from further consideration. For example, there are over 21M gene pairs connected at this level in the GO (BP) branch (similar score = 0.447, Supplementary Fig. 2b ). This led to the removal of 230 genes with no commonly associated MICA with other genes beyond the root term. All ontology-based networks (GO:BP, GO:MF, MPO and HPO) were constructed according to the following procedure summarized in Supplementary Fig. 2c : Pairwise similarity scores given by the procedures above resulted in dense weighted networks. We further applied the disparity filter 19 to extract the backbone of the network and discard structural redundant edges (gene pairs with corresponding disparity p-value > 0.05). The disparity filter provides a dynamic cutoff that considers the strength of the similarity scores of a gene in reference with all similarity values of its neighbors. Similar to using a hard cutoff, edges between gene pairs with low similarity scores (e.g., 06) are virtually unaffected. Edges with medium similarity scores (3