Network-based SNP meta-analysis identifies joint and disjoint genetic features across common human diseases

Shared locus vs. shared association analysis

Despite the size difference, the SLN shows greater network heterogeneity (SLN = 1.30, SVN = 1.17) and lower centralization (SLN = 0.175, SVN = 0.205) and density (SLN = 0.014, SVN = 0.021) values than the SVN. Furthermore, the intersection between the SVN and the SLN lacks not only 5% of the nodes but also 10% of the edges of the SVN. These numbers imply that the process of translating LD data into locus information is at least partly inconsistent. Analysis of the structure of the assigned LD blocks showed two error sources in shared locus analysis (illustrated in Figure2). First, variants in two independent LD blocks are assigned the same locus but are not in LD (assignment error I). Thus, shared loci are found that are not reflected in the variant based data. Second, if two SNPs are in strong LD but the individual long range LD patterns of the SNPs diverge (e.g. the LD block of one SNP covers a greater area at the given r²-threshold), an assignment error II occurs. In this case the two SNPs are assigned to different loci (in the example above, this is due to the different sizes of their LD blocks which may contain distinct gene sets) and their LD connection is lost. These observations suggest that (i), the SLN contains loci which overlap between traits but the associated markers are not in strong LD, (ii), there are several traits which are connected to the SLN via a single, potentially misleading link (as not mirrored in the variant-based data), and, (iii), even a LD-based locus assignment is unable to identify all shared associations (n = 25 of unidentified loci, based on the assignment error II). Due to this limited sensitivity and specificity in detecting LD-based correlations between the reported markers on locus scale, we used the smaller but more accurate variant-based SVN for further analyses. Moreover, the risk allele data in the SVN allows for the inclusion of the direction of the effects, i.e. agonistic and antagonistic, on different traits. Diseases that share antagonistic or agonistic, respectively, associated variants are listed in Tables1 and2. Another advantage is that the SVN can be compared to other variant-based approaches assessing the genetic overlap between traits. Recently, a statistic to identify SNPs with effects across phenotypes (the cross-phenotype meta-analysis statistic, CPMA) was proposed by Cotsapas et al.[8]. It compares the distribution of association P-values of a SNP across seven GWAS on distinct autoimmune diseases to the exponential(1)-distribution representing the expected decay rate of association P-values. As in our approach we use pre-filtered associations, this method cannot readily be employed on our data. However, using the data provided by Cotsapas and colleagues on autoimmune loci in the SVN, we retrieved CPMA P-values on 30 SNPs (∼ 17%) corresponding to 28 loci (∼ 30%) in our data. The CPMA classified all SNPs as significantly effective across diseases (P < 0.05, see Additional file3: Table S2). Thus, we were able to validate nearly one third of the loci contained in the SVN by an independent approach, which underlines the suitability of the network.

Topology of the SVN

Its degree distribution attributes the SVN a scale-free network, i.e. it approximates a power-law (P(k) ∼ k^−γ; γ = 1.32; R² = 0.69) (Additional file4: Figure S2). Interestingly, also when considering the two node types separately, disease nodes (γ = 0.97; R² = 0.71) as well as locus nodes (γ = 2.98; R² = 0.93) show scale-free degree distributions (Additional file4: Figure S2). The scale-free property classifies the network (and its two sets of node types, respectively) as structured, i.e. non-random[13]. It has to be considered that the limited size of the SVN leads to inaccuracies in distribution fitting and thus reduces the explanatory value of this observation. However, as clinically related diseases (i.e. diseases which present similar symptoms) should present a higher genetic overlap than unrelated disorders, this finding meets expectations.

The variant-based SVN also shows no artificial character with regards to its topology. Both locus and disease node sets comprise hubs, here defined as nodes with a degree >3, which form the central elements in the network. As in each GWAS multiple markers are associated with a single disease, one would expect hubs to be constituted mostly of disease nodes. In line with that, 74% of the hubs in our network are disease nodes. The remaining 26% are loci hubs (seven gene loci and three intergenic loci). Several of these loci have been previously identified as influencing susceptibility to multiple diseases like the HLA region on chromosome 6[16], a cancer locus at chromosome 8q24[12], and a coronary artery disease locus at chromosome 9p21[11]. Further hub loci are PTPN22, a known player across several autoimmunity disorders[17], and IL23R, which has been shown to direct inflammatory processes[18]. In addition, we observed hubs which have not yet been described as predisposing to a whole group of diseases, such as TNPO3 which appears to predispose to various autoimmune diseases like systemic lupus erythematosus, systemic scleroderma, and rheumatoid arthritis[19–21], or TNFSF15, which shows associations with several inflammatory diseases[22–25]. As expected, in the majority of cases the traits linked to one hub can be assigned to the same disease group and, further, diseases which are not obviously related to other disorders linked to the respective hub are mostly associated with antagonistic signals. For instance, in a four-gene locus at chromosome 17q12 (GSDML/IKZF3/ORMDL3/ZPBP2, see Additional file2: Table S1), four autoimmune diseases are associated with the same risk allele that in turn has opposite effects on asthma[20, 25–27]. Thus, our results indicate that loci associated with several diseases have an effect specific to a certain disease group rather than effects on unrelated diseases, and that, if there is an effect on an unrelated disease, it can often be distinguished by the direction of the effect.

Disease clustering mirrors trait relatedness

To identify shared and branching mechanisms we split the SNP association data into agonistic and antagonistic variants. Since in most cases there is no solid and comprehensive basis of experimental data that would allow for a more sensitive classification, we suggest that the best available indication of distinct effects of a variant on two diseases is the signal itself being different. Therefore, we define a SNP to be agonistic if all disorders are associated with the same risk allele of the SNP. Conversely, we consider a SNP antagonistic if the associated risk alleles differ between diseases. Accordingly, in the analysis of genetic overlaps as a measure of trait similarity only agonistic variants were included.

As similar diseases are more likely to share associations than diseases in distinct classes, we expected the SVN to be organized in a modular fashion. This was confirmed by the decrease of the degree distribution of the topological coefficient with the number of links per node (Additional file4: Figure S2). To retrieve these modules, we applied a hierarchical clustering approach. The SVN contains two node types (loci and traits). As we wanted to directly assess variant-based disease relatedness, we used its disease centric projection for hierarchical correlation clustering. For data normalization, we calculated the Pearson correlation coefficient (PCC) for all pairs of diseases based on their genetic agonistic overlap with all other diseases. The clustering returned 15 disease clusters (Figure3) and six diseases which show no or only weak correlation with any other disease. With the exception of the heterogeneous cluster 5 (hypertriglyceridemia, ovarian neoplasms, lymphoid leukemia, atopic dermatitis), the clusters mostly contain related diseases. However, many clusters also contain traits unrelated to the other phenotypes like schizophrenia in the autoimmune cluster 2. This indicates that clinical disease classifications appear to be reflected on the genetic level in general. Notably, several small clusters contain diseases which are either linked through common environmental risk factors – like smoking for lung neoplasms, peripheral vascular diseases, and chronic obstructive pulmonary disease – or present high frequencies of comorbidity, e.g. type 2 diabetes (T2D) and obesity. To get an insight into the overall extent of reported comorbidities of the diseases within the 15 clusters, we used publicly available resources[28, 29] and literature mining. The within-cluster fraction of disease co-occurrence ranged from 75% to 100% (μ = 95.89%, σ = 8.66%) which provides empirical evidence of the interrelation of diseases clustered together by genetic information. Such clusters containing diseases that present high ratios of comorbidity may be potential artifacts due to “contaminated” disease cohorts including a substantial number of comorbid cases. The unbiased search for the relation of a trait marker to a disease phenotype as performed in GWAS does not distinguish between markers for a primary or related secondary (comorbid) disease. Our results suggest that this aspect may have been underestimated in some studies, albeit the presence of independently shared etiological mechanisms can naturally not be ruled out in general.

Overall, the outcome of the clustering poses the question of the extent of the influence of phenotype classification and population stratification on GWAS results. Frequent comorbidities (also of seemingly unrelated diseases such as obesity and cancer[30]), diagnostic difficulties in highly related diseases like Crohn’s disease (CD) and ulcerative colitis (UC)[31], and structural (genetic) differences in population subgroups are known to complicate GWAS and impact their outcomes[32]. With growing sample sizes in case–control studies, the potential of false positives produced by such phenomena also increases. As a response, manifold control procedures to handle these and other confounding factors have been developed which are widely used and well appreciated[32]. The heterogeneity of the clusters we retrieved once more highlights the need for the development and application of such methods.

Odds ratio as potential indicator of primary effects

In the context of agonistic association overlap between related diseases, we used the odds ratios (ORs) reported with the SNPs to investigate their impact on the respective traits. In general, the highest ORs are reported for associations of autoimmune diseases to the HLA locus on chromosome 6. Associations with traits where few gene variants with strong effects are reported, e.g. rs6107516 in the prion protein PRNP associated with Creutzfeldt-Jakob disease (OR = 38.5) or rs2071348 in the hemoglobin gene cluster at 11p15.4 associated with beta-thalassemia (OR = 4.33), are exceptions from the majority of associations displaying small ORs[33, 34]. Based on the effect size of variants associated with more than one trait and the same risk allele, we identified three patterns.

First, we identified variants which are likely to present general agonistic risk factors for a group of related diseases or syndromes such as rs13015714 in an interleukine receptor gene cluster at chromosome 2q12.1. This SNP is associated with celiac disease (CeD) and Crohn’s disease (CD) with equal ORs (OR = 1.19)[24, 35]. In cases of frequent comorbidities, though, comparable ORs have limited informative value. The SNP rs9939609 in the FTO gene for instance is associated with T2D and obesity with nearly equal ORs (OR = 1.34 and OR = 1.32) and thus appears to link two coequal traits of the metabolic syndrome (we refer to the definition of the International Diabetes Foundation, 2006)[36, 37]. However, for SNPs in the FTO gene it has been shown that adjustment for body mass index results in the loss of significance (OR∼ 1.0) of the association with T2D[38].

Second, SNPs appeared in several cases to be primarily associated with one disease, which in turn represents a risk factor for another associated trait. For instance, rs2200733 on chromosome 4q25 is linked to atrial fibrillation with a higher OR (OR = 1.72) than to stroke (OR = 1.26)[39, 40]. Another example is rs964184 which is located proximal to the apolipoprotein gene cluster on chromosome 11q23 which is associated with hypertriglyceridemia with a markedly higher OR (OR = 3.28) than to coronary disease (OR = 1.13)[41, 42]. The lower effects of the markers on the hypothesized “secondary sequels” may be explained by the fact that these are caused by the primary diseases, but with less than 100% penetrance.

Third, we speculate that the OR might allow conclusions with respect to the evaluation of an association in cases where similar traits are linked to the same SNP with diverging effect sizes. For instance, CD and ulcerative colitis (UC) share multiple risk loci. The two diseases are strongly related in their etiology and pathology. Thus, a clinical distinction of both diseases is difficult if based only on few criteria and might lead to inaccurate case ascertainment leading to mixed associations[31]. However, for several SNPs such as rs11209026, which is located in the IL23R gene, we found notably higher effects on CD (OR = 3.84) than on UC (OR = 1.74)[25, 43]. Conversely, rs3024505 which lies proximal to the IL10 gene shows a greater effect on UC (OR = 1.46) than on CD (OR = 1.12)[24, 44]. Interestingly, it has been shown that IL23 is selectively upregulated in CD while levels in UC patients are normal and IL10 expression appears to be higher in UC as compared to CD[24, 45]. Thus, the OR might – similar to the above examples – allow for identifying potentially misleading associations in closely related diseases which may result from diagnostic errors.

Identification of branching etiologies

We searched for evidence that antagonistic signals represent genetic indicators of branching points in the etiologies of two diseases or disease groups. For the assessment of potentially multifunctional variants we therefore focused on markers with inverse effects. We identified 44 such variants, which represent almost 4% of the original association data analyzed and about 25% of the SNPs associated with more than one disease. Of those 44 variants, about one fifth (n = 9) are located in the HLA region. SNP-markers in that region are known to differ in their ability to capture the classical HLA-alleles[46] and therefore were not considered further for the present analysis.

For cases where the function of the harboring genes is known, we were able to identify conclusive models. For instance, rs2736100 in the telomerase reverse transcriptase (TERT) gene was reported to exert antagonistic effects in idiopatic pulmonary fibrosis (IPF) and testicular germ cell tumor (TGCT) and two other cancer traits[47–53]. Whereas telomerase activity is generally upregulated in tumors sustaining proliferation and potentiating mutagenesis and transformation of cancer cells[54], in IPF limited cell division due to decreased telomerase activity is thought to contribute to the phenomenon of high percentages of apoptotic cells in fibroblasts[55]. Consistent with that observation, disturbed telomerase activity in TGCT is believed to form a distinct mechanism of cancerogenesis in this tumor type[53]. This distinction from other cancer traits is believed to be based on the fact that testicular germ cells are the only adult cell type with high telomerase expression[56]. Another example is the telomerase RNA component TERC, which is essential for TERT functioning. Opposite alleles of SNP rs10936599 are associated with CeD and colorectal cancer (CRC)[35, 57]. Jones et al. showed that rs2293607, a variant tagged by rs10936599, alone is sufficient to modulate TERC expression[58]. While in CRC this leads to TERC overexpression and longer telomeres, the opposite might apply to CeD, which exhibits telomere reduction and genomic instability[58, 59]. The observation that both constituents of the telomerase complex contain independent antagonistic variants is an intriguing finding. It suggests parallel, autonomous evolution of two functionally interacting loci gone to fixation at a trade-off between early cell senescence or increased apoptosis rates (as in IPF and CeD) and oncogenesis.

A further example is rs1393350 in the tyrosinase (TYR) gene where the opposite alleles are linked to vitiligo and melanoma[60, 61], potentially mirroring the inverse correlation observed for the two traits. The phenomenon is based on the presentation of TYR (self-) antigens on the cell surface of melanocytes. It is hypothesized that in vitiligo the immune system is hypersensitive towards TYR antigens, which are overexpressed in melanoma cells[62]. A possible explanation may be that opposite alleles differentially influence the antigenicity of the TYR protein, thereby conferring protection from melanoma but susceptibility to vitiligo through immune surveillance and vice versa.

In cases of functionally less or uncharacterized genes and their involvement in the associated diseases, our approach can still be used to suggest potential trait-specific effects. Antagonistic effects of rs12720356 (localized in the TYK2 gene) in CD and psoriasis, for instance, might point towards different patterns of cytokine signaling in these two diseases[24, 63, 64]. Likewise, rs12727642 and rs35675666, both located in the PARK7 gene and inversely associated with CeD and UC, could indicate differential effects of oxidative stress on each trait[25, 35].

Variant-based analysis of joint and disjoint genetic features

In this study, we identified overlapping genetic associations and their corresponding loci with analogous or contrasting effects on different diseases. We addressed the methodological challenges of the identification of the functional entities affected by GWAS-detected variants.

Associations formally implicate genomic regions which are captured via tagging SNPs representing haplotype blocks. By using the population-specific LD-based haplotype data provided by the HapMap project[65, 66] or, more recently, the 1000 genomes project[67], SNP arrays are constructed aiming at a high coverage of the total genome variation, but without considering biologically functional aspects. The advantage of GWAS as a method is its unbiased approach to identify genomic regions compromised in a disease; a major drawback is that the association of markers without knowledge of the causal variants and their effects does not allow for a straightforward biological interpretation.

As we show, the reliability of an automated assignment of LD-based loci to the trait-associated variants is strongly context-dependent. Especially in cases of high gene density or, conversely, in intergenic regions/gene deserts, resolving GWAS signals is not possible without further knowledge. Simplifications such as more basic locus assignment approaches which neglect the LD structure of the genome (e.g. classifying a SNP as affecting only the most proximal gene) may seem more intuitive, might facilitate analyses and could be useful to identify causal disease-gene associations. These correct associations of genes which are detected through significant enrichment of a harbored tagging variant in a patient cohort may not be discovered when incorporating LD data in cases where the LD block of the respective variant spans across several genes. However, such approaches disregard a basic principle defining the current GWAS paradigm, namely the use of LD information in the design of genotyping arrays to achieve the genome-wide coverage of common SNPs. Hence, it can be problematic to project the variant-based GWAS data on genes or loci. Accordingly, we decided to use variant-based methods and concentrated on strong gene candidates identified via the gene function of single-gene loci whenever suggesting potential biological effects of the considered variants.

In the analysis of genetic overlaps we followed the hypothesis that the effects of variants shared across several diseases correspond to the reported risk alleles. If the risk allele is the same in all associated diseases, we assume the effect to be the same, i.e. that there is a common underlying etiology. For closely related diseases a positive correlation is not surprising, e.g. a GWAS on psoriatic arthritis (PSA) will also detect agonistic variants such as rs33980500 that are also associated with psoriasis (PS)[68, 69]. Indeed, the vast majority of agonistic variants in our data set links groups of related diseases and thus may mark interesting target regions for therapeutic interventions. However, we also found a few agonistic signals connecting apparently unrelated diseases, e.g. rs6010620 which exerts susceptibility for both glioma and atopic dermatitis (AD)[50, 51, 70, 71]. If our hypothesis is correct, an endophenotype influencing both diseases may be present which has yet to be identified. For antagonistic SNPs, on the other hand, we describe plausible mechanisms that may render variants protective against one trait and predisposing to another, labeling the affected genes/loci as pleiotropic. If pleiotropic effects are as frequent as evolutionary modelers postulate[2, 72] and this effects can be identified by analyses based on GWAS, this might have great implications for the development and use of therapeutics because it would enable avoidance of potential side effects when targeting such loci. Already, there are more than 50 genotype/drug interactions known for which therapeutic dosing recommendations are available[73].

Association selection and curation

We obtained the core list of candidate sentinel SNPs from ‘A Catalog of Published Genome-Wide Association Studies’[1] accessed on June 30, 2011 (http://www.genome.gov/gwastudies). Additional associations where retrieved from HuGE Navigator[76] and automated Text Mining[77]. New (i.e. not contained in the GWAS Catalog) association markers were manually tested on compliance with the criteria for inclusion in the GWAS Catalog before insertion in the candidate list. Associations with copy number variants (CNVs) as well as with pending SNPs were removed. For consistency, SNP identifiers were mapped to the RefSNP numbers of the same dbSNP release (build 131,http://www.ncbi.nlm.nih.gov/snp). For the same reason we semi-automatically translated trait descriptions to the official terms given in the Medical Subject Headings (MeSH). In this process, associations with quantitative and non-disease traits were eliminated. Finally, we lowered the association P-value threshold from 10⁻⁵ (as in the GWAS Catalog criteria) to 10⁻⁷ as to reduce the potential of artificial associations. The workflow of the methods of our approach is sketched in Figure4.

Construction of GWAS networks

For the construction of the locus-based data representation, we defined an associated locus as the whole genomic region captured by SNPs in strong LD, r² ≥ 0.8, with the marker originally reported in a GWAS contained in our data set. The locus is then characterized as all genes located within this genomic region (referred to as “gene locus”) (Figure2). If the region contains no genes, the locus is assigned to its chromosomal location (referred to as “intergenic locus”). LD data and gene information were obtained with the SNAP tool[78]. After locus assignment, our final data set consisted of 111 different traits linked via 1,120 SNPs to 508 gene loci and 226 intergenic loci.

Based on this list we constructed a bipartite graph consisting of two disjoint sets of nodes (Figure1A) representing the complete association data. The first node set corresponds to the traits, whereas the other set comprises the associated loci. Two nodes are connected by an edge if a variant within the respective locus is associated with the corresponding trait. By removal of isolated traits, i.e. traits which share no associated locus with another trait (n = 27) (Figure1A), and cutting out loci which are associated with only one trait (n = 577), we retrieved the SLN (Figure1B).

To obtain a variant-based representation of the data, we repeated the network generation on marker scale by utilizing the set of variants associated with more than one distinct trait. For this, we used the LD data to mutually assign the associated traits of sentinel SNPs in pairwise LD if not already present. In other words, each variant is, in addition to its own associated traits, assigned the traits associated with all correlated SNPs. This set consists of 175 SNPs located in 94 loci and associated with 55 diseases (Additional file2: Table S1). In the resulting bipartite SVN, a trait and a locus are linked if the locus contains a variant which comprises associations with this and at least one other trait. Here, the allele information was included in the graph visualization by coloring of the edges (Figure1C). Both the SLN and the SVN are provided as machine-readable files, see Additional files5 and6.

Network analysis

where N _i is the number of neighbors of i and S(i, j) is the number of shared neighbors of nodes i and j (undefined if i and j do not share a neighbor) plus one if j is a neighbor of i.

for the linear transformation of the power-law functions, i.e. ln y = a + b ln x.

Determination of agonistic and antagonistic effects

For all variants associated with more than one trait, we manually extracted the risk alleles (OR > 1, independently of major or minor allele status) and odds ratios from the reporting studies. The alleles were mapped to the forward DNA strand according to dbSNP 131. The same procedure was applied to markers which were indirectly associated with a trait over LD. If for all traits the same associated risk allele (and corresponding allele, respectively) was reported, the SNP was classified as agonistic. If the risk alleles of a SNP were opposed in the associated diseases, the variant was classified as antagonistic.

Genetic clustering

We applied complete-linkage hierarchical clustering to identify groups of traits genetically overlapping with respect to agonistic signals. Normalization was performed using the linear PCC defined as${\rho }_{X,Y}=\frac{cov\left(X,Y\right)}{{\sigma }_X{\sigma }_Y}$ where the input are the vectors of the variant-based agonistic overlap of two distinct diseases X and Y to all other diseases. Thus, disorders which are clustered together show a homogeneous association overlap pattern to all other diseases, while diseases which are not clearly assigned to a cluster present a more heterogeneous pattern relatively unique in the SNP data. For cluster definition, we used a Euclidian distance threshold of 1.71. This threshold was determined as the maximal distance at which the six traits not correlating with other diseases (Figure3) remain non-clustered.

Calculation of the CPMA statistic for autoimmune loci

We downloaded the dataset S1 from[8] and extracted the information on autoimmune-linked SNPs contained in the SVN. We used the Z-scores given in the file to compute two-sided P-values for all seven GWAS. Using the CPMA code provided onhttp://www.cotsapaslab.info/index.php/software/cpma/ we calculated the CPMA P-values as described in[8].