A large-scale analysis of tissue-specific pathology and gene expression of human disease genes and complexes

Systematic Generation of an Atlas of Disease-Associated Protein Complexes with Tissue Resolution.

By mining the GeneCards (15) resource for genes associated with diseases, we generated a list of 2,227 unique disease-related proteins. Similar to the method that we reported earlier (13), an in silico approach for generating disease-associated protein complexes based on an inferred human protein–protein interaction network was used [see supporting information (SI) Text and Fig. S1]. Following this strategy, we generated 1,524 raw complexes comprising 45,662 unique interactions between 5,202 unique proteins. The quality of the complexes was validated by measures identical to the ones reported in major experimental screens in Saccharomyces cerevisiae, Escherhichia coli, and Homo sapiens (7–11, 16, 17), showing that the quality of our data matches the reproducibility, average probabilistic interaction scores, accuracy, and coverage reported in these studies and that the complexes are true biological entities (see SI Text, Table S1, and Fig. S2). Finally, the complexes were mapped to tissues by using the expression data from 73 nondiseased tissues from the Novartis Research Foundation Gene Expression Database (GNF) (18). The expression level of a complex in a tissue was calculated by averaging over the expression levels of all genes represented in the complex.

Mapping Complexes to Diseases.

To map complexes to diseases we systematically identified the proteins that had been associated to each of the diseases mentioned in OMIM. This was done by using the protein to OMIM mapping displayed in GeneCards (http://www-bimas.cit.nih.gov/cards/) database. We then measured the overlap between proteins in complexes and proteins associated with the diseases and calculated the significance of this overlap. Because a number of complexes are known to be involved in different diseases we allowed for a complex to be associated with more than one disease. In total the 1,524 raw complexes were mapped to 1,054 OMIM diseases. In the further text we refer to these as disease complexes.

Disease–Tissue Association Matrix.

To our knowledge no systematic mapping of diseases to affected tissues exists. We determined the covariance of a disease with a tissue by identifying the number of publications comentioning the disease and tissue (and synonyms thereof), relative to the number of publications mentioning the disease or tissue alone (19). We transformed the covariance into an association score between a tissue and a disease by calculating the fraction of covariance that a given tissue–disease pair constituted, of the total covariance for a given disease. Calculating an association score for the 73 tissues used in the GNF tissue atlas (18) versus 1,054 OMIM diseases yielded a disease–tissue association matrix (Fig. 2). By manually validating the associations we determined a cutoff where tissues associated with the pathology of a given disease could be determined with a precision of >80% (see SI Text, Table S2, and Fig. S3), meaning that above this threshold tissues relevant to the pathology of a given disease can be accurately identified among the GNF atlas tissues in >80% of the cases. High confidence associations scoring above this threshold are blue to dark blue in Fig. 2. Tissues associated with the pathology of a given diseases are in the further text defined as disease–tissue associations scoring above this cutoff.

Mapping Complexes to Cancers.

A large number of genes have been associated with cancers, due to aberrant expression or somatic mutations in tumors. However, few of these genes have actually been proven to play a role in the initiation of the tumor. Hence, an automated mapping of cancer genes to complexes would include many genes that are mutated in tumors but do not cause the cancer. Because we are interested in studying the tissue distribution of disease-initiating genes and complexes, we manually created an exhaustive list of heritable cancer genes that initiate tumors through germ-line mutations. These genes were mapped to OMIM diseases describing the cancers manually (Table S3). For this subset of genes, there is compelling evidence that defects are the primary cause of the cancer. In total we extracted a subset of 51 genes in which mutations lead to heritable cancers and mapped them to 59 cancers. Because most cancer mutations are either loss or gain of function that could influence the mechanisms of disease progression and have bearing on the mechanisms of tissue specificity, we further stratified the cancer genes into loss or gain of function as defined in Vogelstein et al. (4). Examples of loss-of-function genes are tumor suppressor or DNA repair genes that become defective when mutated, and examples of gain of function are kinases that become constitutively activated by mutations (Table S4). Cancer-associated complexes were identified as complexes enriched for this subset of genes. In the further text we refer to these as cancer complexes.

Generating a Disease–Tissue Association Matrix for Cancers.

Cancer to tissue association mapping is not straightforward. In this study we were interested in exclusively studying the tissues in which tumors are initiated through germ-line mutations of particular genes. Because cancers generally affect many tissues through downstream effects such as metastases, associations to noninitiating tissues had to be filtered out. Furthermore, many cancer syndromes, arising from germ-line mutations in cancer genes, also include nonmalignant pathology, for which disease–tissue association had to be disregarded in this analysis. For this reason, we manually analyzed the complete subset of tissues associated with heritable cancer syndromes resulting in a precision approximating 100% for the cancer–tissue associations (SI Text and Table S5).

Correlation Between Pathology and Tissue-Specific Expression.

First, we analyzed the expression of disease genes in the tissue with the highest disease association in the disease–tissue matrix (rank 1). This analysis was repeated for the 2nd to 25th highest associated tissues (rank 2 to 25) and the average z score at each rank level was plotted as a curve (Fig. 3A). For example, myosin heavy chain 6 (MYH6) is involved in hypertrophic cardiomyopathy and the tissues from the GNF atlas ranked first and second in relation to hypertrophic cardiomyopathy are heart and cardiac myocytes. We determined the z score of MYH6 in heart (tissue rank 1), the average z score of MYH6 in the 2 highest ranked tissues, heart and cardiac myocytes (tissue rank 2). This procedure is repeated for ranks 3 to 25. This gives a set of rank-dependent z scores for MYH6. This procedure is repeated for every disease gene in every disease yielding rank-dependent z scores for every gene–disease combination, which is plotted in Fig. 3A. This figure shows the clear tendency of overexpression for disease genes in tissues with the highest rank (blue curve). The curves for cancer genes show 2 different trends. Although gain-of-function genes are overexpressed in tissues with the highest rank (red curve), loss-of-function genes are underexpressed (green curve).

To see whether the observed expression trends were significant, we averaged the z scores in the tissues associated with the disease and compared the scores with their expression levels in nonaffected tissues (Fig. 3B). For non-cancer disease genes we observed a significant tendency of overexpression (P < 1.0E-6), which is also the case for gain-of-function cancer genes (P = 3.9E-2), but with less significance. Loss-of-function cancer genes show the converse trend of underexpression (P = 1.0E-2).

We carried out the same analysis for the protein complexes which showed that the expression trend observed for disease genes is conserved at the level of disease protein complexes (see Fig. 3 D and C). These disease complexes display a significant tendency to be overexpressed in tissues where they are involved in pathology (P < 10E-6, blue curve). While protein complexes significantly enriched for gain-of-function cancer genes follow the tendency of overexpression (P = 0.44, red curve), complexes enriched for loss-of-function cancer genes are underexpressed (P = 3.4E-3, green curve).

Because the z scores were lower for the cancer genes and complexes compared with the more robust values of the non-cancer disease genes and complexes, we tested whether this result was influenced by the dataset and normalization method. We replicated the analysis by using a different robust multiarray (RMA)-based normalization scheme (20). Expression data normalized with this algorithm still showed a significant overexpression of disease genes and complexes, but both the over and underexpression trends for the cancer genes and complexes decreased in significance. To test whether a few diseases or tissues were driving the observed trend, we analyzed the expression trend broken down into single tissues (Fig. S4a) and by bootstrapping the dataset both on disease and tissue level. This analysis shows that most tissues and diseases contribute to the observed results and they are robust to bootstrapping of the dataset (Fig. S4b).

Examples of Disease Complexes with Tissue and Phenotype Correlation.

Examples of the correlations found between tissue expression and pathology or phenotype reported are provided in Fig. 4. Also, the most significant gene ontology (GO) subcelluar and functional categories for the complex in question are indicated followed by the significance with which the complex can be assigned to this GO category (21). Tissue names are as defined in the GNF atlas. The full sets of proteins in each complex can be seen in Fig. S5.

XY sex reversal can be caused by mutations in the transcription factors SRY (Sex determining Region Y) (22), SOX 9 (the SRY sex determining region Y-box 9 gene) (23), NR5A1 (the nuclear receptor subfamily 5A1), more commonly known as SF1 (24, 25); and NR0B1 (nuclear receptor subfamily 0B1), more commonly known as DAX1 (26). Additionally, SOX 9 is associated with campomelic dysplasia, a bone disorder that leads to a number of associated skeletal and cartilaginous deformities (27). SF1 is needed for gonad and adrenal differentiation (25, 28) and for proper steroidogenesis as well as for Mullerian Inhibiting Substance (MIS) ligand and MIS receptor expression (28, 29). DAX1, which leads to XY sex reversal both when overexpressed, by inhibiting SF1 (26), and when inactivated, as it is required for testis differentiation by regulating expression of SOX9 (30). Although the activity of SF1, DAX1, and SOX9 is required for testis differentiation, development, and maintenance, none of these genes are essential for ovarian development and maintenence (30–33). Here, we identify a transcriptional regulation complex (GO:0006355: P = 1.9E-8) containing DAX1, SF1, and SOX9, all of which are known to be associated with sex reversal (P = 6.9E-6). Furthermore, the complex contains SOX8 that is closely related to SOX9 and implicated in regulating the expression of testis-specific genes (34). Whereas the complex is overexpressed in testis cells, it is underexpressed in ovaries (Fig. 4), which coincides with the known biology of the most well characterized of its components. Our method thus has predictive value because it can (i) detect interactions between molecules that, by themselves, are known to be important in sex differentiation and determination by producing sex reversal, (ii) validate these findings by demonstrating dimorphic tissue-specific expression that correlates with the pathology resulting from inactivation of several members of the complex, and (iii) reveal the importance of new interactors worthy of further study.

Four other complexes, where tissue-specific overexpression correlates with pathological manifestations, are depicted in Fig. 4 (see SI Text and Fig. S5 for more details on these 4 complexes and for examples of cancer-related complexes). These include (i) a complex involved in Charcot–Marie–Tooth disease type 4F and overexpressed in spinal cord, dorsal root ganglion, and skeletal muscles; (ii) a sarcoglycan complex involved in Limb–girdle muscular dystropy overexpressed in skeletal muscle, cardiac myocytes, and heart; (iii) a myofibril complex involved in familial cardiomyopathy overexpressed in several tissues associated with the disease such as heart and cardiac myocytes; (iv) and a complex involved in catechol metabolism and Parkinson disease, overexpressed in a number of relevant brain tissues including the caudate nucleus, subthalamic nucleus, and globus pallidus. Although the overexpression of the sarcoglycan and myofibril complex in muscle tissues is well known, the ovarian–testes dimorphic expression pattern of the sex-reversal complex, and the overexpression of a Parkinson complex in several relevant brain tissues of the basal ganglia are suggestive of the underlying tissue-specific biology of these disorders. Across all examples the tissue-specific expression patterns correlate with the pathological changes observed when one or several members of the complex are defective.

Mapping Genes and Complexes to Tissues.

We used the GNF tissue atlas (18) that includes reproduced RNA expression experiments from 79 human tissues. Six tissues were removed because they were derived from cancer tissues. We chose the GNF dataset because it displays high reproducibility (44), and the transcript levels show generally a linear relationship with protein abundance (39). We log-transformed hybridization levels and normalized within each tissue (to ensure equal weight), followed by a normalization across all tissues, thereby ensuring that expression levels represented the relative presence of a transcript in one tissue compared with the other 72 healthy tissues in the dataset. For complexes, the normalized expression levels of all genes in a complex were averaged for each tissue. To test the effect of different normalization methods on our results, we prepared the same dataset with Eklund and Szallasi's normalization method (20) and compared the results.

A Curated Set of Genes in Which Mutations Lead to Tumor Formation.

We curated a set of genes in which mutations had been shown to lead to heritable tumor formation and mapped them to OMIM diseases (see Table S3). By following the definitions introduced by Vogelstein et al. (4) we also noted whether the genes were oncogenes or nononcogenes (such as tumor suppressors or DNA repair proteins) (see Table S4).

Mapping of Complexes to OMIM Diseases.

We calculated the enrichment of proteins involved in the same OMIM disease by using the annotations in GeneCards, which has previously been shown to be an accurate way of mapping genes to diseases (13). We calculated the significance of an enrichment by using a hypergeometric test.

Under- and Overexpression Significance.

We averaged the expression z score over all disease genes in the most disease-associated tissue as determined from the disease–tissue matrix. For each rank from 1 through 25, we calculated the average z score yielding a curve. In Fig. 3A, this curve is plotted as the average z scores of all gene–disease pairs in tissues with a particular rank. This procedure was repeated for gain-of-function and loss-of-function cancer genes. Again this approach was repeated on a protein complex level. All reported significances are 2-tailed using the Student's t test.

Disease-Tissue Association Matrix.

To identify the tissues most affected by diseases described in the OMIM database (14), we used comentioning of a given disease with a given tissue across PubMed (19). The tissue names from the Novartis Research Foundation Gene Expression Database (GNF) (18) were manually curated and translated to corresponding medical subject heading (MeSH) terms (to reduce errors in word recognition and handle synonyms properly). Similarly, the disease names were determined by using disease titles provided in OMIM. Also, these titles were manually curated and translated to the relevant MeSH terms. We used Ochiai's coefficient (OC) as a measure of similarity derived from the cooccurrences (45–47), and calculated an association score (see below), as the percentage of the total normalized cooccurrence of a given disease that could be attributed to a given tissue. Validation was carried out as described in SI Text.