In this study, we develop and apply a novel algorithm that generates tissue-specific functional networks in the laboratory mouse by integrating diverse functional genomic data, and we demonstrate that our tissue-specific networks are more accurate in predicting phenotype-related genes than a single global functional network. In the following sections, we first outline the strategy used to generate tissue-specific networks by interrogating gene expression profiles across tissues and integrating different data sources using Bayesian statistics. Second, we developed a cross-network comparison metric for identifying significantly changed genes across networks which are enriched in tissue specification and development. Third, we quantitatively demonstrate that combining our tissue-specific networks with a state-of-the-art machine learning algorithm can produce improved predictions of genotype-phenotype relationships compared to previous single global networks 
. Fourth, we identify candidate genes related to male fertility specifically predicted by our tissue-specific networks (but not by the global network), and verify a top prediction in an independent, unbiased mutant screen. Finally, we used cerebellum-specific network to predict genes associated to a less-studied disease, ataxia, which are supported by both literature and experimental evidence. The predictions made by our approach for all examined networks are available online at http://mouseMAP.princeton.edu
Constructing a Bayesian model to establish tissue-specific functional networks
A common mechanism resulting in tissue-specific protein functionality is the modulation of gene expression levels between tissues 
. This observation is our theoretical foundation for establishing tissue-specific networks, in which links between proteins represent the probability that they are involved in the same biological processes within a specific tissue. To simulate such tissue-specificity, we developed a Bayesian approach () that incorporates highly-reliable, low-throughput measures of tissue-specific gene expression into training set, which we utilized to produce networks focused on the real functional relationships occurring within the tissue under consideration. This Bayesian framework essentially learns how informative each dataset is given a set of ‘gold standard’ training pairs, i.e.
pairs of proteins known to be functional in the same biological process and both expressed in the tissue of interest.
Strategy for constructing tissue-specific networks and predicting phenotype-associated genes.
In the global (non-tissue-specific) sense, following previous definitions 
, ‘gold standard positives’ are defined by co-annotation to specific Gene Ontology (GO) biological process terms 
, while ‘gold standard negatives’ are defined as pairs that both have specific GO annotations yet do not share any annotations. For each tissue-specific gold standard set, a positive pair has to meet two requirements: first, the pair must be ‘co-functional’ as defined in the global sense, and second, both genes must be expressed in the tissue under consideration as evident in highly reliable, low-throughput expression datasets, which, in most cases, is necessary for the pair to have a functional relationship in that tissue. These tissue-specific gold standards are then used to quantify how relevant each genomic dataset is in recovering tissue-specific functional relationships, regardless of the tissue of origin for each genomic dataset. This allows us to leverage the entire compendium of high-throughput genomic data to generate accurate tissue-specific networks, even for tissues which do not have existing tissue-specific whole-genome experiments, by relying on non-tissue-specific datasets, heterogeneous samples, and potentially related tissues and experiments. For example, biliary tract, which is not specifically represented in our current collection of high-throughput features used for classification, can still be accurately predicted by utilizing information from related, heterogeneous samples, such as gene expression microarrays of whole liver or the hepatic system, as well as non-tissue-specific information, such as sequence phylogeny and in vitro
binding assays. Thus our approach can leverage the implicit relationships between a dataset and a tissue and therefore enables generation of tissue-specific networks even from feature data that is not resolved for a specific tissue type.
For tissue-specific expression information, our gold standards rely on the Gene Expression Database (GXD) of the Mouse Genome Informatics group (MGI). GXD provides an extensive, hierarchically structured dictionary of anatomical expression results for mouse to allow us to carry out our analysis 
. The data in GXD are derived from traditional, “small-scale” expression experiments, such as in situ
hybridization, RT-PCR, and immunohistochemistry, which simply reflect presence or absence of a gene within the tissue examined. No high-throughput expression data were used for our gold standard construction. In total, there are 107 tissues included in our analysis.
We pursue two main goals in this study: First, we generate tissue-specific networks that synthesize as much data as possible and provide these networks to the public through an online visualization interface at http://mouseMAP.princeton.edu
. For this, we gathered diverse genomic data for mouse as inputs (Dataset S1
) to support the functional relationships, including protein-protein physical interactions 
, homologous functional relationship predictions from simpler organisms 
, phenotype and disease data 
and 960 expression datasets, totaling 13632 experimental conditions 
. The reliability of each dataset is learned through Bayesian network classifier training, using the tissue-specific gold standards described above. Essentially, a dataset deemed more relevant and accurate for the tissue under consideration will be given higher weight, and the final probability of pair-wise functional relationships is determined by updating the initial probability (prior) based on the weighted input of all genomic datasets. This procedure resulted in tissue-specific probabilistic functional relationship networks for the laboratory mouse that effectively summarize these diverse data sources and enable biology researchers to easily explore the resulting functional landscape. Second, we test the hypothesis that tissue-specific networks could assist us to predict phenotype-related genes more accurately. In this case, to prevent circularity in our methodology, phenotype and disease data were excluded from network generation, and the results were used to predict novel phenotype-associated candidate genes. We demonstrate that tissue-specific networks enhance biological clarity and result in more accurate predictions. Our resulting networks and predictions provide biology researchers with functional interactions specific to each tissue as well as phenotype hypotheses of genes.
Robust recovery of tissue-specific functional relationships
One key application of tissue-specific networks is to identify novel genes and relationships between genes that may be specific to a particular tissue. To computationally evaluate our ability to identify novel relationships, we used cross-validation to test whether our tissue-specific Bayesian scheme is more accurate than the global network. Cross-validation was used to assess predictions by evaluating the accuracy of recovering subsets of known annotations withheld during the training process. Specifically, we performed 3-fold cross-validation, by holding out one third of the tissue-specific ‘gold standard’ pairs in each of the three iterations. We learned the parameters in the Bayesian networks, i.e., the reliability of each dataset, through the other two thirds of the ‘gold standard’, and then used these networks to predict the probabilities for the held-out one third of the protein pairs.
Compared to a single global functional relationship network, our approach significantly improved our ability to predict tissue-specific functional linkages. The mean AUC (area under the receiver operating characteristic curve, which represents the accuracy in recovering tissue-specific functional relationships) for the global network estimated through three-fold cross-validation was 0.68. Tissue-specific networks achieved median AUC of 0.72. With a random baseline of 0.5 in AUC, this represents a ~20% improvement of the tissue-specific networks over the predictive power of the global network. This improvement is consistent over all 12 major organ systems defined by GXD 
. (). Immune system-related networks acquired the most median improvement of 22.7% and digestive system-related networks achieved least median improvement of 14.3%. For example, for lymphoid system (MA:0002435), we improved our AUC from 0.65 to 0.72 and for ventricular zone, brain, we improved from 0.65 to 0.77. Such improvement is consistent across the entire precision-recall spaces (, Dataset S2
for all precision-recall curves). In all cases, tissue-specific networks performed better than the global network in predicting functional relationships specific to that tissue, which demonstrates the robustness of our integration approach across different systems and tissues in the laboratory mouse.
Tissue-specific networks are more accurate than the global network in reflecting protein functional relationships.
Network comparison reveals altered gene-gene functional connections across tissues enriched for activated biological processes
One important application of our tissue-specific networks is to identify functional relationships between genes that change significantly across tissues. This provides a platform for analyzing tissue-specific molecular interactions, as well as tissue-specific roles for genes that are ubiquitously expressed but play different roles in different tissues. For example, Wnt10b (wingless related MMTV integration site 10b) is expressed in many tissues throughout development and participates in many biological processes including bone trabecular formation 
and cell differentiation involved in skeletal muscle development 
. The interactors of Wnt10b in our muscle-specific and bone-specific functional networks reflect its differential roles in these two tissues. The top neighbors in the muscle-specific network consist of genes responsible for skeletal muscle development (). For example, BIN1 participate in the biological process muscle cell differentiation (GO:0042692) 
, PLAU is involved in the process skeletal muscle tissue regeneration (GO:0043403) in rat and MYF6 directly function in muscle cell biogenesis 
. In fact, 8 out of the 19 top connected nodes of Wnt10b in the muscle-specific network are involved in skeletal muscle cell development, reflecting the functional role of Wnt6b. On the contrary, in the bone-specific functional network, the top neighbors of Wnt10b consist of genes involved in bone mineralization and bone structure formation (), representing 12 out of the 19 top connected nodes. This observation suggests that our networks can provide a resource for comparing the dynamic functions of a single gene across different tissues.
Top connected genes of Wnt10b in muscle-specific and bone-specific networks.
To quantify gene connectivity changes across networks, we developed a metric that captures how much the edges involving a gene differ across networks (see methods
), and we implemented a web-based visualization interface (http://mouseMAP.princeton.edu
) allowing users to query genes of interest and compare the local network between tissues. Essentially, connectivity change of a gene is defined by the sum of absolute values of fold changes (over prior) of connections between this gene to all other genes. Some genes vary greatly in their connectivity between tissues, potentially reflecting their tissue-specific roles. Of the top 100 altered genes, they were significantly enriched for “anatomical structure development” (GO:0048856) and “organ development” (GO:0048513). Additionally, genes with connectivity altered in specific tissues compared to the global network, tend to be enriched for GO terms related to the tissue under consideration. For example, when comparing the nervous system-specific network (MA:0000016) against the global network, the top changed genes are enriched in “central nervous system development” (GO:0007417), “diencephalon development” (GO:0021536), and “brain development” (GO:0007420) (). The full enrichment analysis is provided in Dataset S3
Example enriched Gene Ontology terms in the tissue MA:0000016 nervous system.
Prioritizing phenotype-associated genes in relevant tissue-specific networks
A key hypothesis in this study is that analyzing tissue-specific networks may improve our ability to identify phenotype-related genes. To test this hypothesis, we regenerated tissue-specific networks using the same Bayesian approach as above, but excluded all phenotype and disease data as inputs to avoid circularity in our cross-validations. Then, we mapped 451 phenotypes to their most related tissue in the laboratory mouse according to the terminology and description of these phenotypes in the Mammalian Phenotype ontology 
. For each phenotype, we compared novel predictions made using the appropriate tissue-specific network as compared to using the global network. This method is based on our previously developed machine learning scheme (network-based SVM) 
that mines information in functional relationship networks to prioritize candidate genes according to their links to known genes related to a disease or phenotype.
To test whether our tissue-specific networks are more capable of identifying phenotype-associated genes than the global network, we used bootstrap bagging 
to evaluate which network performs better. Bootstrap bagging is suitable for phenotype predictions, where positive examples (known phenotype-associated genes) and negative examples (random genes) are highly imbalanced 
. Its stability and comparably good performance in estimating error rates has been tested in extensive simulations for positive example set sizes ranging from less than 20 
to >200 
, which is the approximate range we are using in our evaluation. For the 451 mapped phenotypes, the median AUC when utilizing tissue-specific networks is 0.794, representing an improvement of 11.8% over utilizing the global functional network. For many phenotypes, using tissue-specific networks can improve our ability to extract potentially experimentally-verifiable predictions. For example, at one percent recall (the low recall end is where most of the follow-up experimental confirmations will focus on), we achieved a precision of 1.00 compared to 0.33 using global network for the phenotype abnormal spleen white pulp morphology (MP:0002357), and a precision of 0.5 compared to 0.28 for abnormal malpighian tuft morphology (MP:0005325). Additionally, the AUC for “abnormal osteogenesis” (MP:0000057) was 0.77 using the global network, but 0.81 for tissue-specific networks. The AUC for “abnormal nervous system electrophysiology” (MP:0002272) using the global network was 0.716, but was 0.763 using the nervous system-specific network ( for example precision-recall curves). Such significant improvement demonstrates the potential of mining tissue-specific networks to prioritize phenotype-associated genes.
Tissue-specific networks perform better than the global network in predicting genes related to different phenotypes.
Performance improvements were consistent across phenotypes of different sizes (). For phenotypes with 300–1000 annotated genes (around 1.5% to 5% of genome), we achieved a median AUC of 0.814 (improvement of 8.7%); for phenotypes with 100–300 genes, the median AUC was 0.792 (improvement of 13.0%); and for phenotypes with 30–100 genes, the median AUC was 0.769 (improvement of 11.0%). At 10 percent recall for the 300–1000, 100–300, and 30–100 groups, we achieved precisions of 14.8, 17, and 20 fold over random, respectively. This consistency indicates the robustness of tissue-specific networks against the number of known genes in predicting phenotype-associated genes.
Performance improvements were also consistent across different major organ systems. Phenotypes involved in the endo/exocrine system achieved the most significant improvement in AUC (+35%, compared to global networks against baseline of 0.5) and those in cardiovascular system achieved 21.8% improvement in AUC. However, prediction accuracy was improved across all major systems, with the least improvement of 5.9% in renal/urinary phenotypes. Phenotypes related to musculoskeletal systems achieved the highest AUC of 0.82 and the group with lowest AUC was digestive system, which still achieved an average of 0.78. The consistency in improvements across different organ systems demonstrates the robustness of our modeling framework to predict phenotype-related genes in a tissue-specific manner.
Predicting and testing phenotype/disease genes using the tissue-specific networks
We focused on two cases to illustrate how our tissue-specific networks can facilitate disease gene discovery. These two phenotypes represent two extremes of the phenotype/disease-associated gene prediction problem. The first, reduced male fertility, is a broadly defined, common phenotype with many causative genes already known. The second, ataxia, is a rare neurological disorder affecting ~3–10/100,000 of the general population 
. Roughly 40 genes are known to be associated with this disease, but the majority of both familial and sporadic cases remain unexplained. Predicting candidate genes related to rare genetic diseases is challenging in that little prior knowledge is available for these diseases. These phenotypes are related to two different tissue-categories (reproductive and neurological systems), enabling us to highlight the broad applications of our approach across organ systems. We used these two examples and experimental confirmations to demonstrate the power of tissue-specific networks to discover disease genes.
First, we used male fertility related phenotypes to test the performance of tissue-specific networks to predict phenotype-related genes. To do so, we utilized a recent, nearly comprehensive literature review of genes involved in mammalian spermatogenesis and male fertility phenotypes 
, which we organized into a hierarchy of male fertility-related phenotypes (Dataset S4
). This curation effort is independent of, and more comprehensive than, the current GO or MP annotations related to male fertility, which makes these lists excellent, non-circular test sets. We tested whether the testis-specific network could predict male fertility genes more robustly than the globally integrated network, and found that the testis-specific network significantly improved our ability to predict spermatogenesis-related phenotypes. For example, for predicting genes related to ‘spermatid head and nuclear modifications,’ we achieved 4.5-fold improvement in precision at 1 percent recall; for ‘acrosome-related genes,’ we achieved 3.6-fold improvement; and for ‘germ/Sertoli cell interaction genes,’ we achieved 3.3-fold improvement. On the other hand, for terms that are not specifically related to male-reproductive systems, such as ‘association with methylation and acetylation,’ and ‘association with Golgi Apparatus,’ we observe no performance improvements using the testis-specific network. This illustrates that tissue-specific functional relationship networks are tuned to predict phenotypes closely related to these tissues.
We selected Mybl1
to demonstrate the specific utility of the male-reproductive network to predict fertility related genes. Mybl1
(MGI:99925) is among our top candidates in multiple phenotypes related to male fertility, including ‘association with chromatoid body and manchette’, ‘transcription factor involved in spermatogenesis’ and, ‘spermatogenesis’. However, in the global network, Mybl1
was not a strong candidate for these phenotypes, as it was predicted with negative values. Therefore Mybl1
is an ideal candidate to test the accuracy of our tissue-specific network-based phenotype predictions. In our male-reproductive network, the majority of the top interactors of Mybl1
are indeed well-known male fertility genes (), including Dmc1
(required for meiosis and male fertility 
(a DEAD-box helicase required for male, but not female, germ cell development 
(encoding testis-specific cytochrome c 
) and Lhx9
(a LIM homeobox required for sex differentiation and normal fertility 
). Moreover, Mybl1
was independently identified recently in an unbiased mutagenesis screen for infertility phenotypes involving meiotic arrest 
. We found that the Mybl1
mutants are characterized by low testis weight and depletion of male germ cells, as shown in . Additionally, analysis of the mutant testis transcriptome suggested that MYBL1 is a “master regulator” of the meiotic cell cycle and transcriptional program 
, and at least one gene regulated by MYBL1, Cyct
, is among the top interactors of MYBL1 predicted by our network. Together, these findings on an infertility phenotype and suggestions of a corresponding potential mechanism confirm the accuracy of predictions from our tissue-specific network and show that when taken with expression analyses and other data, they can be used as a basis for functional testing.
Prediction and verification of infertility-related genes through male reproductive system-specific networks.
In addition to the well-studied phenotype of male infertility, we also examined a less well-understood disease, ataxia, to investigate whether our tissue-specific networks can identify genes related to phenotypes or diseases with limited prior knowledge. Gene identification through genetic approaches, such as pedigree analyses, has had a major impact on our understanding of ataxia (over 40 candidate genes identified so far). Genetic testing is now an integral part of assessment. Routinely, a blood sample of any new ataxia case is mailed in for laboratory evaluation. However, the majority of the sporadic cases as well as the familial cases are so far unexplained. We curated the known gene list (43 in total) related to human ataxia, mapped these genes to their mouse orthologs, and used this list as seeds to predict additional candidate genes using our cerebellum-specific network, which is the major tissue affected by ataxia.
Our cerebellum-specific network reveals connections of ataxia-related genes not shown in the global network. A key, known ataxia gene is Atcay
(ataxia, cerebellar, Cayman type homolog (human)), and in the cerebellum-specific network, two of its top interactors are Cacna1e
(with connection confidence 0.943, ranked 18) and Grm1
(0.902, ranked 46) (). These are plausible candidate genes since Grm1
is a known mouse ataxia gene 
, and Cacna1e
encodes a subunit of an R-type calcium channel, while mutations to the related protein family member Cacna1a
, encoding a subunit of an L-type calcium channel, causing spinocerebellar ataxia. However, in the global network these interactions are much weaker (0.647 for Cacna1e
and 0.763 for Grm1
respectively), and would not be identified in the top 100 connections of Atcay
, which supports the utility of tissue-specific networks relevant to ataxia to identify candidate genes ().
Top connected genes to Atcay in the cerebellum-specific network reveals likely ataxia candidates.
In addition to identifying these novel, likely correct edges, we also identified novel candidates using our SVM-based approach described above. Out of our top 10 novel candidates, we found strong evidence in the literature for 4 of these genes to be associated to ataxia (); suggesting at least a 40% success rate at low levels of recall. Among these, SORBS1 physically interacts with ATXN7, an autosomal dominant gene causing cerebellar ataxia 
. RBFOX1 physically interacts with the c-terminus of ATXN2, another autosomal dominant gene causing cerebellar ataxia 
. It is thought that RBFOX1 might contribute to the restricted pathology of spinocerebellar ataxia type 2 (SCA2) 
. The homozygous mouse knockout of a third gene, Plcb4
induces ataxia, although no human patients have been identified with mutations in this gene. A fourth gene, Plp1
, is implicated in Spastic paraplegia-2 and Pelizaeus-Merzbacher diseases 
, which are disorders closely related to ataxia. It is also a homologue of Pmp22
, which is involved in Charcot Marie Tooth disease type 1A, a sensory neuropathy common in some forms of ataxia 
. Thus, even in the case of less well-studied phenotypes or diseases, our tissue-specific approach is able to identify likely candidates as evidenced by our success rate of at least 40% for ataxia-related predictions based on the cerebellar network, compared to a background detection rate of less than 1/500.
Evidence for top 10 predictions for ataxia-causing genes using mouse cerebellum-specific networks.