Skip to main content

iLOCi: a SNP interaction prioritization technique for detecting epistasis in genome-wide association studies

Abstract

Background

Genome-wide association studies (GWAS) do not provide a full account of the heritability of genetic diseases since gene-gene interactions, also known as epistasis are not considered in single locus GWAS. To address this problem, a considerable number of methods have been developed for identifying disease-associated gene-gene interactions. However, these methods typically fail to identify interacting markers explaining more of the disease heritability over single locus GWAS, since many of the interactions significant for disease are obscured by uninformative marker interactions e.g., linkage disequilibrium (LD).

Results

In this study, we present a novel SNP interaction prioritization algorithm, named iLOCi (Interacting Loci). This algorithm accounts for marker dependencies separately in case and control groups. Disease-associated interactions are then prioritized according to a novel ranking score calculated from the difference in marker dependencies for every possible pair between case and control groups. The analysis of a typical GWAS dataset can be completed in less than a day on a standard workstation with parallel processing capability. The proposed framework was validated using simulated data and applied to real GWAS datasets using the Wellcome Trust Case Control Consortium (WTCCC) data. The results from simulated data showed the ability of iLOCi to identify various types of gene-gene interactions, especially for high-order interaction. From the WTCCC data, we found that among the top ranked interacting SNP pairs, several mapped to genes previously known to be associated with disease, and interestingly, other previously unreported genes with biologically related roles.

Conclusion

iLOCi is a powerful tool for uncovering true disease interacting markers and thus can provide a more complete understanding of the genetic basis underlying complex disease. The program is available for download at http://www4a.biotec.or.th/GI/tools/iloci.

Background

A major challenge for human genetics is identifying susceptibility genes for complex heritable diseases. Advanced single nucleotide polymorphism (SNP) genotyping technology and genome-wide association study (GWAS) are at the forefront of research in this area. In conventional single locus analysis, each variant is tested individually for disease association. Systematic analysis of GWAS data in this manner can typically uncover multiple SNPs associated with complex diseases [13]. These analyses have provided valuable insights into the genetics of complex diseases; however, they typically detect only common, low-risk variants each with small effect and explain only a tiny proportion of disease heritability [4].

The existence of interactions among genes (epistasis) has been proposed to constitute a major proportion of disease heritability, which is not captured by single-locus GWAS [5]. The genetical nature of epistasis can be described by several different models as shown in a variety of interaction schema discussed in [6]. Note that genetic factors primarily function through a complex mechanism; thus, epistatic interactions are not limited to independent gene pairs. Multiple genes interacting through a biological network (i.e. indirect interactions) exist which can modify disease penetrance and expressivity.

A number of methods for detecting epistatic interactions among genotypic data have been proposed. Most methods employ a statistical approach to identify interacting marker pairs based on deviation from a null distribution and estimation of type I error. These statistical approaches have been shown to work well in theory, e.g., regression methods [7, 8], partitioning chi-square [9], Focused Interaction Testing Framework (FITF) [10], Bayesian model selection [11], and other recent approaches [12, 13]. However, the need for control of type I error reduces power to detect interactions in real data, which is exacerbated by the huge number of statistical tests performed in this analysis [14].

Given the challenges for statistical approaches, non-statistical methods such as machine-learning and data-mining methods have been proposed for the study of genetic interactions [15, 16]. Instead of model fitting, these methods attempt to explain all of the heritability in terms of marker interactions. Multifactor dimensionality reduction (MDR) is an brute-force method for identifying the most plausible interactions which fit the data [17]. However, MDR and other recently published exhaustive non-parametric approaches [18] are computationally complex and thus impractical for analysis of GWAS data. To overcome the computational burden of non-parametric analysis, several techniques have been developed that employ statistics to assist the non-parametric search for epistasis, including SNPHarvester [19], SNPRuler [20], and BOOST [21]. In these methods, the search space is reduced by a filtering step, usually employing a statistical threshold. The filtered dataset is then used for non-parametric search for epistasis. Although these methods can be applied for analysis of GWAS data, the interactions found rarely offer any new insights since the majority of interacting markers map to the same genomic regions. For example, the analysis of WTCCC (Wellcome Trust Case Control Consortium) data by BOOST revealed that after removal of linked pairs, no interactions were found for five of the seven diseases. Using another approach for exhaustive search of interactions, the most recent paper by Ueki and Tamiya [22] also reported very few interactions in the WTCCC data.

The possible reason for the disappointingly modest improvement of the current hybrid approaches is that they do not adequately account for marker dependencies not related to disease. A well known marker dependency which can confound the identification of genomic regions associated with disease is linkage disequilibrium (LD). LD is non-random association of genotypes at two or more loci that can be on the same or different chromosomes. LD is caused by a number of factors, including genetic linkage and the rate of recombination [23]. Earlier reports [24, 25] showed that LD contrast, i.e., differences in LD patterns between case and control groups can reveal the disease signal above the noise of background LD in candidate disease regions. However, to our knowledge, LD contrast has not been employed for comprehensive genetic epistasis study, owing to the high computational complexity.

Clearly, a computationally efficient and comprehensive prioritization technique is required which accounts for marker dependencies unrelated to disease. Moreover, instead of trying to control type I error, a prioritization procedure may be more effective in revealing more of the true disease markers which may have modest individual effects and interact in complex higher-order networks.

In this paper, we propose a novel tool for prioritizing gene-gene interactions called iLOCi (interacting Loci). The iLOCi algorithm ranks all SNP pair combinations according to a novel heuristic that we call ρdiff. The iLOCi program is specifically designed to handle large-scale GWAS data partly through the application of data parallelization. The tests with WTCCC datasets show that the top ranked pairs by our algorithm reveal novel disease genes, several of which are consistent with biological networks underpining disease etiology.

Methods

iLOCi algorithm

The proposed iLOCi algorithm performs genome-wide analysis for identifying SNP pairs that are plausibly associated with a disease. No prior genetical assumptions are employed in the algorithm, which allows the exploration of different dimensions of the association results. The framework can be characterized into two main modules: 1) calculating SNP pair dependencies separately in case and control groups and 2) disease SNP pair prioritization as shown in Figure 1.

Figure 1
figure 1

The workflow of the iLOCi algorithm.

Calculation of SNP pair dependencies

iLOCi explores all possible combinations of SNP pairs. Given N SNPs from a SNP array with the SNP index starting from 1 to N, there are a total of ( N 2 ) = N ( N 1 ) 2 possible pairs. Each SNP pair is assigned a unique index (i,j), where i≠j.

From the large number of SNP pairs, it is necessary to identify the dependency unrelated to disease. This dependency includes linkage disequilibrium (LD), population structure, genotype calling artifacts, etc. and is performed separately between the case and control groups. This step of the algorithm is called dependence test. Therefore, for each indexed SNP pair, the algorithm calculates two scores, ρ case and ρ control . The calculated ρ values using genotypic information were proven to be concordant with LD values (see Additional file 1). LD values are calculated using allelic deviation from the Hardy-Weinberg Equilibrium (HWE) model, which assumes that, without the introduction of specific disturbing factors, the frequencies of alleles and genotypes in a population remain constant from one generation to the next. However, it should be noted that the only information captured by ρ values is the correlation between markers, which is needed for identifying interactions. For LD calculation, the haplotypic phase is also considered, which is computationally very demanding for datasets of this size.

To compute marker ρ values, each SNP locus is considered as a discrete random variable and the numeric values of -1, 0 and 1 are assigned to homozygous wild (w), heterozygous (h), and homozygous variant (v) types respectively. This encoding ensures zero-means, which obviates a normalization step. Let x and y be two discrete random variables of SNPx and SNPy, respectively. Let P(x,y)represents a genotypic joint probability mass function, whose entries are the probability of genotype combinations from both SNPs. Hence, there are nine possible genotypic combinations that are represented by the following matrix:

P ( x , y ) = [ P w w P w h P w v P h w P h h P h v P v w P v h P v v ]

For example, P ww is a probability that (x,y) are both homozygous wild type. Each of these probabilities can be calculated by dividing the number of the joint genotypic outcomes with the total number of individuals for either case (Ncase) or control (Ncontrol) groups. For example, P w w ctrl = P ( x = w , y = w ) ctrl = N ( x = w , y = w ) c t r l N ctrl . The dependence test must be performed for all possible SNP pairs. The correlation value ρcontrol for each SNP pair is calculated as:

= { x w y w P w w ctrl + x w y h P w h ctrl + x w y v P w v ctrl } + { x h y w P h w ctrl + x h y h P h h ctrl + x h y v P h v ctrl } + { x v y w P v w ctrl + x v y h P v h ctrl + x v y v P v v ctrl } ( x w 2 P x = w ctrl + x h 2 P x = h ctrl + x v 2 P x = v ctrl ) × ( y w 2 P y = w ctrl + y h 2 P y = h ctrl + y v 2 P y = v ctrl ) = { ( 1 ) ( 1 ) P w w ctrl + ( 1 ) ( 0 ) P w h ctrl + ( 1 ) ( 1 ) P w v ctrl } + { ( 0 ) ( 1 ) P h w ctrl + ( 0 ) ( 0 ) P h h ctrl + ( 0 ) ( 1 ) P h v ctrl } + { ( 1 ) ( 1 ) P v w ctrl + ( 1 ) ( 0 ) P v h ctrl + ( 1 ) ( 1 ) P v v ctrl } ( ( 1 ) 2 P x = w ctrl + ( 0 ) 2 P x = h ctrl + ( 1 ) 2 P x = v ctrl ) × ( ( 1 ) 2 P y = w ctrl + ( 0 ) 2 P y = h ctrl + ( 1 ) 2 P y = v ctrl ) = P w w ctrl P w v ctrl P v w ctrl + P v v ctrl ( P x = w ctrl + P x = v ctrl ) × ( P y = w ctrl + P y = v ctrl )

Note that P x = w ctrl , P x = v ctrl , P y = w ctrl , P y = v ctrl are the estimated probability of SNPx wild type, SNPx variant type, SNPy wild type and SNPy variant type respectively.

By the same reasoning, ρcase is calculated as:

= P w w case - P w v case - P v w case + P v v case ( P x = w case + P x = v case ) × ( P y = w case + P y = v case )

Disease SNP pair prioritization

The next step is to identify whether the same SNP pair (x,y) from case and control groups have contrasting patterns of ρ values. A difference test is performed by differentiating the ρ values between the case and control groups using a simple subtraction operation, namely ρdiff=|ρcontrol-ρcase|.

To select the highly associated SNP pairs, all SNP pairs are ranked according to the ρdiff values. The ranking of top SNP pairs was chosen, rather than a P- value cutoff in order to avoid too many false positive pairs due to the heavy-tailed distribution phenomenon, where the Gaussian distribution decreases faster than the distribution of disease associated SNP pairs [26].

Parallel computing algorithm implemented in iLOCi

The iLOCi algorithm is designed for genome-scale analysis which requires the computation of a huge number of SNP interaction pairs, e.g.≈1.25x1011 pairs for a 500,000 SNP dataset. Data parallelization is applied to accelerate this computationally intensive and time-consuming process. The SNP interaction matrix is divided into submatrices of 100,000 or fewer SNPs each. Each SNP interaction submatrix is computed in parallel using a MacPro workstation with 2×2.4 GHz quad-core Intel Xeon processors with 8GB RAM. With this configuration, the complete WTCCC dataset can be analyzed in 19 hours. Details for implemention of the code and data parallelization are available upon request.

Testing iLOCi algorithm performance using simulated data

The performance of iLOCi for detecting disease-associated gene interactions was evaluated and compared with FastEpistasis [27]. The evaluation was made using simulated datasets, which were generated using the GenomeSIM program [28]. The algorithm performance was determined for detection of four different epistatic interaction scenarios:

  1. 1)

    Single pair interaction without marginal effects: Eighteen epistatic models in [29] with heritability (h2) of 0.2, 0.3, and 0.4 were used for performance comparison (see Additional file 2: Table S1). These heritability levels were chosen to represent those typically found in common complex diseases. The minor allele frequency (MAF), which is the frequency of the less common allele, was assigned to be two levels, 0.2 and 0.4. In total, there are six model groups comprising three models with the same heritability and MAF for each group. 100 independent datasets containing 1600 samples (800 cases and 800 controls) with 100 SNPs were generated for each model group.

  1. 2)

    Single pair interaction with marginal effects: Six epistatic models in [30] with MAF of 0.5 were tested (see Additional file 2: Table S2). 100 independent datasets containing 800 samples (400 cases and 400 controls) and 100 independent datasets containing 1600 samples (800 cases and 800 controls) with 100 SNPs each were generated for each model group.

  2. 3)

    Multiple independent interacting pairs without marginal effects: Eight models of multiple interactions described in supplementary material of [19] were tested. Each of these models were generated from five epistatic models described in [29]. Each model used the same heritability and MAF. 100 independent datasets containing 1600 samples (800 cases and 800 controls) and 100 SNPs were generated for each model group.

  3. 4)

    Higher-order interactions: Data were simulated for the eight interaction network models based on pairwise interaction described in [31] for three-, four-, and five-loci interating networks (see Additional file 2: Table S3). 100 independent datasets containing 800 samples (400 cases and 400 controls) were generated. The number of SNPs varies from model to model.

The algorithm performance was demonstrated by the percentage of accuracy, which is determined by the proportion of 100 independent datasets in which the algorithm correctly identified the interacting SNP pairs. For situations 1 and 2, the identification of disease SNP pair is defined as correct if the disease SNP pair is the top ranked pair with the highest ρdiff score (for iLOCi) or the lowest P-value (for FastEpistasis). For multiple independent interacting pairs (case 3), the identification is taken as correct when all five disease SNPs fall in the top five ranked pairs with highest ρdiff score (for iLOCi) or lowest P-value (for FastEpistasis). The prediction of higher-order interactions is defined as correct when all disease SNPs are found within all top ranked pairs. The top ranked pairs are defined as all consecutive pairs comprising at least one disease SNP in each pair.

Testing algorithm performance using the WTCCC dataset

In addition to the simulated data, our algorithm was applied to the real genotypic data of WTCCC (Wellcome Trust Case Control Consortium) [3]. This dataset encompasses ~500,000 SNP genotypic data of ~17,000 British samples which are divided into 3000 shared control samples and ~2000 case samples for each of seven complex diseases: bipolar disorder (BD), coronary artery disease (CAD), Crohn's disease (CD), hypertension (HT), rheumatoid arthritis (RA), type1 (T1D) and type2 (T2D) diabetes.

For these real datasets, data cleaning was required prior to the analysis. We considered only SNPs and individuals passing WTCCC data quality control [3]. We further filtered the SNP set using MAF>0.05 leaving 355,882 SNPs (complete set) for all diseases. We also generated a SNP marker gene-only subset of 176,148 present in genes (defined as within 10Kb flanking an annotated gene model reported in RefSeq version 36.3).

First, ρdiff values for the seven WTCCC diseases were calculated for all possible (≈63×109 for complete and ≈15×109 for the gene-only subset) pairs. Next, the empirical ρdiff distributions for each disease were graphed using kernel density plot. For the gene-only SNP subset analysis, the top ranked 1000 SNP pairs were chosen for functional analysis to uncover biological significance. From these pairs, a list of genes was extracted based upon RefSeq (version 36.3) physical locations of SNPs in the genome. To understand the biological significance of the novel genes reported by our algorithm, we also used the candidate gene prioritization feature of ToppGene [32] using the cutoff of P-value = 0.01 with Bonferroni correction. The training sets for the ToppGene candidate gene prioritization were the lists of all genes reported in the HuGE Navigator database [33] for the seven diseases. The test sets for the ToppGene analysis were the lists of novel (not reported in HuGE Navigator database) genes represented among the top ranked 1000 SNP pairs obtained from iLOCi.

Results

iLOCi algorithm validation

We used simulated datasets to validate the iLOCi algorithm for identifying various disease-associated epistatic interactions. We chose FastEpistasis for performance comparison with iLOCi due to the fact that the data were simulated according to an interaction model; hence this tool would be most suitable for testing. Moreover, the theoretical basis for FastEpistasis is widely accepted for genome-wide analysis.

The first result testing for a single interacting pair demonstrated that the top ranked iLOCi pair was the disease interacting pair in 18 different inheritance models without the presence of marginal effects. Overall, its performance was approximately the same as FastEpistasis for most of the model groups and slightly better in some cases (h2=0.2, MAF = 0.4; h2=0.3, MAF = 0.4) as shown in Figure 2. For epistatic interactions with marginal effects, iLOCi outperformed FastEpistasis in most models, except in model 2 for which both methods failed to detect the interacting disease marker pair (Figure 3). Furthermore, we want to demonstrate the specificity as well as sensitivity of iLOCi for detecting multiple interacting disease marker pairs as would be present in a real dataset. Therefore, the receiver operating characteristics (ROC) were plotted for different thresholds of ranked marker pairs, and for different models of heritability and MAF (Figure 4). Generally, iLOCi has high sensitivity and specificity, although the performance tends to be worse with lower degrees of heritability. Moreover, it should be noted that the minimum ρdiff scores that give 100% sensitivity vary greatly from 0.00511 to 0.41663.

Figure 2
figure 2

The performance comparison between iLOCi (I) and FastEpistasis (F) on epistatic models without marginal effects. The algorithm performance is shown as the percentage of accuracy, which is the number of simulated datasets (out of 100) in which the correct SNP pair is identified. The accuracy was tested for two different MAF (0.2, 0.4) and three different levels of heritability (A) 0.4, (B) 0.3, and (C) 0.2.

Figure 3
figure 3

The performance comparison between iLOCi and FastEpistasis on epistatic models with marginal effects. The percentage of accuracy is shown for two different sample sizes (800 and 1600) for six different pairwise interaction models (A-F).

Figure 4
figure 4

The Receiver Operating Characteristic (ROC) curves for simulation datasets of hybrid models. The ROC curves are displayed for five independent interacting SNP pairs.The MAF and heritability parameters were varied (A) h2=0.1, MAF = 0.2, (B) h2=0.2, MAF = 0.2, (C) h2=0.3, MAF = 0.2, (D) h2=0.4, MAF = 0.2, (E) h2=0.2, MAF = 0.4, (F) h2=0.1, MAF = 0.4. The ρdiffvalues are shown that give the maximum true positive rate with the lowest false positive rate (red dots).

In addition to independent interacting pairs, we examined the ability of iLOCi and FastEpistasis to detect higher-order interactions of 3, 4, and 5 loci disease interaction networks for eight models at each level (Figure 5). iLOCi can detect all eight models for all levels of interactions; however, FastEpistasis failed to identify all S3 model interactions. Furthermore, FastEpistasis could detect, with higher than 50% accuracy, in fewer than 50% of the 4-loci network models and only Ep1, Ep3 and Ep5 of the 5-loci network models.

Figure 5
figure 5

The performance comparisons between iLOCi and FastEpistasis on high-order interaction models. The percentage of accuracy is shown for different models (Ep1, Ep3, Ep5, Ep6, Het1, Het3, S1, S3) of high-order interactions among (A) three-loci, (B) four-loci, and (C) five-loci.

In conclusion, these experiments with simulated data validated the iLOCi algorithm for identifying all four types of higher-order gene interaction. iLOCi performance was comparable to FastEpistasis for a variety of two-locus interaction models; however, iLOCi was markedly superior for detecting high-order interactions. This would be a major advantage of iLOCi for analysis of real data since high-order interaction is the type of interaction likely to be found in real data of complex diseases and may account for current missing heritability.

iLOCi analyses of WTCCC data

The iLOCi algorithm was tested against real data obtained from WTCCC. The distribution of ρdiff values follows a Weibull distribution pattern for all seven diseases (Figure 6). From the Weibull distribution with k = 1 and λ=0.018, we calculated P-values for ρdiff scores ranging from 0.05 to 1.0 (see Table 1). For the seven diseases, we selected the top 1000 pairs for which the calculated minimum P-values vary from <2.22e-16 to 1.14e-7 in complete SNP set analysis, and from <2.22e-16 to 4.72e-5 in gene-only SNP analysis (see Table 2).

Figure 6
figure 6

The frequency distribution of ρ diff values from WTCCC datasets. The plot shows the empirical probabilty density function (pdf) generated from combined ρdiffvalues from all seven diseases of WTCCC datasets.The pdf plots generated from each disease are indistinguishable from combined pdf. The plots for Weibull distribution (k = 1, λ=0.018) and Chi-square distribution (degree of freedom = 1) are shown in the same axes.

Table 1 The lookup table of P-values for the associated ρdiff scores
Table 2 The ρdiff scores of the 1st and 1000th ranked SNP pairs and their associated P-values

From iLOCi analysis using the complete SNP marker set, it was found that the great majority of the SNPs have not been previously reported to be associated with the diseases [3]. Furthermore, the majority of these SNPs also do not map to annotated genes. The list of top 1000 SNP pairs is available in Additional File 3. For each disease, iLOCi identified 'hub' SNPs, i.e. SNPs that pair with many other SNPs, e.g., rs1553460 pairs with 1000 other SNPs in BD (Table 3).

Table 3 The hub SNPs/genes identified in the top-ranked 1000 SNP pairs

Owing to the fact that the majority of interacting SNPs do not map to annotated genes, we re-analyzed the data using the gene-only SNP subset. 'Hub' SNPs were also observed at the gene level (Table 3). From this analysis, it was noted that the top ranked 1000 SNP pairs of all seven diseases map to 321 disease-gene associations that have been annotated on the HuGE Navigator database (see Table 4, Additional File 4). On the other hand, the majority of the disease interacting genes among these pairs reported by iLOCi are novel. Moreover, most of these genes were not reported in the original WTCCC study (Table 4). To evaluate the biological significance of the novel genes among these pairs, the ToppGene candidate gene prioritization tool was employed. The full results are shown in Additional Files 3 and 4. Among the novel genes identified by iLOCi, it was observed that some well known disease pathways from KEGG [34] contain several of these genes (see Additional File 5). For instance, the 'neuroactive ligand-receptor interaction' pathway in BD contains 4 novel genes in addition to 11 previously reported genes (Figure 7). Other prominent disease pathways include 'cytokine-cytokine receptor interaction' for CAD (Figure 8) and 'type I diabetes mellitus' for T1D (Figure 9).

Table 4 The disease association of iLOCi selected genes from gene-only SNP analyses
Figure 7
figure 7

The iLOCi detected genes from BD and their maps in 'neuroactive ligand-receptor interaction' pathway. The KEGG pathway diagram [34] shows the mapping of BD-associated genes identified among 1000 top ranked iLOCi pairs in 'neuroactive ligand-receptor interaction' KEGG pathway. The gene families containing the genes previously reported in HuGE Navigator database and the novel disease genes are highlighted in the red boxes and the blue boxes, respectively, with their associated gene names.

Figure 8
figure 8

The iLOCi detected genes from CAD and their maps in 'cytokine-cytokine receptor interaction' pathway. The KEGG pathway diagram [34] shows the mapping of CAD-associated genes identified among 1000 top ranked iLOCi pairs in 'cytokine-cytokine receptor interaction' pathway. The gene families containing the genes previously reported in the HuGE Navigator database and novel disease genes are highlighted in the red boxes and the blue boxes, respectively, with their associated gene names.

Figure 9
figure 9

The iLOCi detected genes from T1D and their maps in 'type I diabetes mellitus' pathway. The KEGG pathway diagram [34] shows the mapping of T1D-associated genes identified among 1000 top ranked iLOCi pairs in the 'type I diabetes mellitus' KEGG pathway. The gene families containing the genes previously reported in HuGE Navigator database and the novel disease genes are highlighted in the red boxes and the blue boxes, respectively, with their associated gene names.

Discussion

In this study, we have developed a new pairwise SNP-interaction prioritization algorithm for GWAS. We hypothesized that by first accounting for pairwise marker dependencies among case and control groups, it would be possible to observe true disease interactions above the noise of dependent markers unrelated to disease, as was proposed in earlier studies of LD contrast (see Background).

In GWAS data, it is well known that LD generates strong pairwise dependency signals that are used to identify disease associated SNPs by imputation. However, this type of signal predominates pairwise markers in analysis of gene interactions. For example, in the approach used by Wan et al. [21], the majority of the interactions identified for all seven WTCCC datasets can be attributed to LD effect, i.e., the interacting SNPs are within 1Mb of each other in the same genomic region. To validate our approach correcting for pairwise dependencies unrelated to disease SNP interactions, extensive tests were performed on simulated data. For a simple model with only one interacting pair, the top ranked iLOCi pair is correctly identified as the disease marker pair. When testing for multiple interacting pairs, iLOCi has high accuracy under the conditions of high heritability and informativeness, i.e., low MAF. On the other hand, low heritability and/or informativeness leads to type I error as observed by ROC plot. In general, the ρdiff scores reflect the degree of heritability and informativeness. Hence, it is not possible to use a single ρdiff cutoff for identifying disease interactions in the real case when the heritability and informativeness are unknown.

From analyses of real GWAS data, it was found that the ρdiff distributions for all seven diseases could be represented by a single kernel density function with Weibull distribution. However, the range of ρdiff values varies among the diseases and follow the known heritability pattern, i.e., HT has the lowest heritability and lowest top ρdiff score, while T1D has the highest heritability and highest top ρdiff score (Table 2). Although it is possible to calculate P-values of the interacting pairs and use them as cutoffs for prioritization, we consider the use of P-value cutoffs inappropriate. For example, a P-value of 1e-5 (corresponding to ρdiff values of approximately 0.2 or greater) would give approximately 16 million significant pairs for T1D and 200,000 pairs for HT. The same phenomenon of unacceptable type I error was found by others when using FastEpistasis for analysis of real datasets. It is debatable whether Bonferroni correction is valid since the tests are not independent, as shown by the heavy-tailed distributions of ρdiff . Current methods for correction of type I error by false discovery rate are also likely to be impractical because of the requirement for permutation testing.

Instead of using P-value significance thresholds, we used the top ranked 1000 SNP pairs for prioritization, which account for a very small portion (<0.0001%) of all possible pairs. Rather than attempting to identify all gene interactions, which practically can not be found [35], we limit the prioritization to the top ranked pairs that are most likely to contain the genetic interactions which are informative of the disease etiology, i.e., disease pathways. From the full SNP set analysis, several hub SNPs were identified for each disease which interact with many other SNPs. For some diseases such as T1D, these hub SNPs map to well-known disease associated genes. However, hub SNPs for BD, HT, and CD do not map to genes. These hub SNPs may mediate interactions at an unknown gene regulatory level, e.g. as non-coding RNAs, miRNAs or cis-regulatory elements. Since our knowledge of gene regulation is far from complete [36], we repeated the iLOCi analysis on the gene-only SNPs subset. By restricting the analysis to SNP pairs in genes only, the ToppGene systems approach for gene prioritization was appropriate, as used by others for GWAS data [3739].

Gene-based prioritization of the interacting SNP pairs revealed significant representation of previously described disease associated genes. Therefore, we are confident that the novel genes found among the prioritized SNP pairs are novel disease-associated genes. For each disease, hub genes were found which pair with many other genes. Some of these disease hub genes are known and have been replicated as disease genes by conventional single-SNP GWAS, including the MHC gene HLADQB1 for T1D and TCF7L2 for T2D. However, some hub genes have not been reported previously, e.g. the CACNG1 gene for RA. This gene's SNP shows a modest P-value (>1e-4) for association by single SNP analysis [3]; therefore, the disease association of this SNP is dependent on multiple interactions with other loci. For each disease, including those with low heritability such as HT, we are able to suggest novel genes and pathways for further investigation, including re-analysis of other GWAS datasets for the same diseases.

Conclusions

In this article, we introduce a novel SNP interaction prioritization method, called iLOCi. The algorithm is computationally efficient, and thus suitable for exhaustive search for interactions along markers in a typical GWAS dataset. We have shown that the approach taken by iLOCi in which marker dependencies unrelated to disease are accounted for reveal genetic interactions of biological relevance to complex disease.

References

  1. Easton DF, Pooley KA, Dunning AM, Pharoah PD, Thompson D, Ballinger DG, Struewing JP, Morrison J, Field H, Luben R, et al: Genome-wide association study identifies novel breast cancer susceptibility loci. Nature. 2007, 447 (7148): 1087-1093. 10.1038/nature05887.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  2. Samani NJ, Erdmann J, Hall AS, Hengstenberg C, Mangino M, Mayer B, Dixon RJ, Meitinger T, Braund P, Wichmann HE, et al: Genomewide association analysis of coronary artery disease. N Engl J Med. 2007, 357 (5): 443-453. 10.1056/NEJMoa072366.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  3. The Wellcome Trust Case Control Consortium: Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature. 2007, 447 (7145): 661-678. 10.1038/nature05911.

    Article  PubMed Central  Google Scholar 

  4. Manolio TA, Brooks LD, Collins FS: A HapMap harvest of insights into the genetics of common disease. J Clin Invest. 2008, 118 (5): 1590-1605. 10.1172/JCI34772.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  5. Moore JH, Asselbergs FW, Williams SM: Bioinformatics challenges for genome-wide association studies. Bioinformatics. 2010, 26 (4): 445-455. 10.1093/bioinformatics/btp713.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  6. Musani SK, Shriner D, Liu N, Feng R, Coffey CS, Yi N, Tiwari HK, Allison DB: Detection of gene × gene interactions in genome-wide association studies of human population data. Hum Hered. 2007, 63 (2): 67-84. 10.1159/000099179.

    Article  CAS  PubMed  Google Scholar 

  7. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, et al: PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007, 81 (3): 559-575. 10.1086/519795.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  8. Zhao J, Jin L, Xiong M: Test for interaction between two unlinked loci. Am J Hum Genet. 2006, 79 (5): 831-845. 10.1086/508571.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  9. Yang Y, Houle AM, Letendre J, Richter A: RET Gly691Ser mutation is associated with primary vesicoureteral reflux in the French-Canadian population from Quebec. Hum Mutat. 2008, 29 (5): 695-702. 10.1002/humu.20705.

    Article  CAS  PubMed  Google Scholar 

  10. Millstein J, Conti DV, Gilliland FD, Gauderman WJ: A testing framework for identifying susceptibility genes in the presence of epistasis. Am J Hum Genet. 2006, 78 (1): 15-27. 10.1086/498850.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  11. Zhang Y, Liu JS: Bayesian inference of epistatic interactions in case-control studies. Nat Genet. 2007, 39 (9): 1167-1173. 10.1038/ng2110.

    Article  CAS  PubMed  Google Scholar 

  12. Ueki M, Cordell HJ: Improved statistics for genome-wide interaction analysis. PLoS Genet. 2012, 8 (4): e1002625-10.1371/journal.pgen.1002625.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  13. Wu X, Dong H, Luo L, Zhu Y, Peng G, Reveille JD, Xiong M: A novel statistic for genome-wide interaction analysis. PLoS Genet. 2010, 6 (9): e1001131-10.1371/journal.pgen.1001131.

    Article  PubMed Central  PubMed  Google Scholar 

  14. Hunter DJ, Kraft P: Drinking from the fire hose--statistical issues in genomewide association studies. N Engl J Med. 2007, 357 (5): 436-439. 10.1056/NEJMp078120.

    Article  CAS  PubMed  Google Scholar 

  15. Cordell HJ: Detecting gene-gene interactions that underlie human diseases. Nat Rev Genet. 2009, 10 (6): 392-404.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  16. McKinney BA, Reif DM, Ritchie MD, Moore JH: Machine learning for detecting gene-gene interactions: a review. Appl Bioinformatics. 2006, 5 (2): 77-88. 10.2165/00822942-200605020-00002.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  17. Ritchie MD, Hahn LW, Roodi N, Bailey LR, Dupont WD, Parl FF, Moore JH: Multifactor-dimensionality reduction reveals high-order interactions among estrogen-metabolism genes in sporadic breast cancer. Am J Hum Genet. 2001, 69 (1): 138-147. 10.1086/321276.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  18. Yoshida M, Koike A: SNPInterForest: a new method for detecting epistatic interactions. BMC Bioinformatics. 2011, 12: 469-10.1186/1471-2105-12-469.

    Article  PubMed Central  PubMed  Google Scholar 

  19. Yang C, He Z, Wan X, Yang Q, Xue H, Yu W: SNPHarvester: a filtering-based approach for detecting epistatic interactions in genome-wide association studies. Bioinformatics. 2009, 25 (4): 504-511. 10.1093/bioinformatics/btn652.

    Article  CAS  PubMed  Google Scholar 

  20. Wan X, Yang C, Yang Q, Xue H, Tang NL, Yu W: Predictive rule inference for epistatic interaction detection in genome-wide association studies. Bioinformatics. 2010, 26 (1): 30-37. 10.1093/bioinformatics/btp622.

    Article  CAS  PubMed  Google Scholar 

  21. Wan X, Yang C, Yang Q, Xue H, Fan X, Tang NL, Yu W: BOOST: A fast approach to detecting gene-gene interactions in genome-wide case-control studies. Am J Hum Genet. 2010, 87 (3): 325-340. 10.1016/j.ajhg.2010.07.021.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  22. Ueki M, Tamiya G: Ultrahigh-dimensional variable selection method for whole-genome gene-gene interaction analysis. BMC Bioinformatics. 2012, 13 (1): 72-10.1186/1471-2105-13-72.

    Article  PubMed Central  PubMed  Google Scholar 

  23. Hedrick PW: Genetics of populations. 2005, Sudbury, Boston, Toronto, London, Singapore: Jones and Bartlett Publishers, 3

    Google Scholar 

  24. Wang T, Zhu X, Elston RC: Improving power in contrasting linkage-disequilibrium patterns between cases and controls. Am J Hum Genet. 2007, 80 (5): 911-920. 10.1086/516794.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  25. Zaykin DV, Meng Z, Ehm MG: Contrasting linkage-disequilibrium patterns between cases and controls as a novel association-mapping method. Am J Hum Genet. 2006, 78 (5): 737-746. 10.1086/503710.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  26. Embrechts P, Klüppelberg C, Mikosch T (eds.): Modelling Extremal Events for Insurance and Finance. 1997, Berlin: Springer Verlag, 1

    Chapter  Google Scholar 

  27. Schupbach T, Xenarios I, Bergmann S, Kapur K: FastEpistasis: a high performance computing solution for quantitative trait epistasis. Bioinformatics. 2010, 26 (11): 1468-1469. 10.1093/bioinformatics/btq147.

    Article  PubMed Central  PubMed  Google Scholar 

  28. Dudek SM, Motsinger AA, Velez DR, Williams SM, Ritchie MD: Data simulation software for whole-genome association and other studies in human genetics. Pac Symp Biocomput. 2006, 499-510.

    Google Scholar 

  29. Velez DR, White BC, Motsinger AA, Bush WS, Ritchie MD, Williams SM, Moore JH: A balanced accuracy function for epistasis modeling in imbalanced datasets using multifactor dimensionality reduction. Genet Epidemiol. 2007, 31 (4): 306-315. 10.1002/gepi.20211.

    Article  PubMed  Google Scholar 

  30. Moore J, Hahn L, Ritchie M, Thornton T, White B: Application of genetic algorithms to the discovery of complex models for simulation studies in human genetics. Proceedings of the Genetic and Evolutionary Computation Conference: July 9-13, 2002 2002; New York, USA. 2002, Morgan Kaufman, 1150-1155.

    Google Scholar 

  31. Neuman RJ, Rice JP: Two-locus models of disease. Genet Epidemiol. 1992, 9: 347-365. 10.1002/gepi.1370090506.

    Article  CAS  PubMed  Google Scholar 

  32. Chen J, Bardes EE, Aronow BJ, Jegga AG: ToppGene Suite for gene list enrichment analysis and candidate gene prioritization. Nucleic Acids Res. 2009, W305-311. 37 Web Server

  33. Yu W, Gwinn M, Clyne M, Yesupriya A, Khoury MJ: A navigator for human genome epidemiology. Nat Genet. 2008, 40 (2): 124-125. 10.1038/ng0208-124.

    Article  CAS  PubMed  Google Scholar 

  34. Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M: KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 2012, D109-114. 40 Database

  35. Zuk O, Hechter E, Sunyaev SR, Lander ES: The mystery of missing heritability: Genetic interactions create phantom heritability. Proc Natl Acad Sci USA. 2012, 109 (4): 1193-1198. 10.1073/pnas.1119675109.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Esteller M: Non-coding RNAs in human disease. Nat Rev Genet. 2011, 12 (12): 861-874. 10.1038/nrg3074.

    Article  CAS  PubMed  Google Scholar 

  37. Dick DM, Aliev F, Krueger RF, Edwards A, Agrawal A, Lynskey M, Lin P, Schuckit M, Hesselbrock V, Nurnberger J, et al: Genome-wide association study of conduct disorder symptomatology. Mol Psychiatry. 2010, 16 (8): 800-808.

    Article  PubMed Central  PubMed  Google Scholar 

  38. Edwards AC, Aliev F, Bierut LJ, Bucholz KK, Edenberg H, Hesselbrock V, Kramer J, Kuperman S, Nurnberger JI, Schuckit MA, et al: Genome-wide association study of comorbid depressive syndrome and alcohol dependence. Psychiatr Genet. 2012, 22 (1): 31-41. 10.1097/YPG.0b013e32834acd07.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  39. Lascorz J, Forsti A, Chen B, Buch S, Steinke V, Rahner N, Holinski-Feder E, Morak M, Schackert HK, Gorgens H, et al: Genome-wide association study for colorectal cancer identifies risk polymorphisms in German familial cases and implicates MAPK signalling pathways in disease susceptibility. Carcinogenesis. 2010, 31 (9): 1612-1619. 10.1093/carcin/bgq146.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

JP is supported by the new researcher grant from the Thailand Research Fund and National Center for Genetic Engineering and Biotechnology (grant number TRG5580011). ST would like to acknowledge the TRF grant number RSA5480026 and the Research Chair Grant 2011 from the National Science and Technology Development Agency (NSTDA), Thailand that partially support this work. ST was supported in part by the office of the higher education commission and Mahidol University under the national research university initiative. This study makes use of data generated by the Wellcome Trust Case-Control Consortium. A full list of the investigators who contributed to the generation of the data is available from http://www.wtccc.org.uk. Funding for the project was provided by the Wellcome Trust under award 076113 and 085475.

This article has been published as part of BMC Genomics Volume 13 Supplement 7, 2012: Eleventh International Conference on Bioinformatics (InCoB2012): Computational Biology. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcgenomics/supplements/13/S7.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Sissades Tongsima.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

JP designed the algorithm and the experiments, generated simulated data, analyzed test results, and wrote the manuscript. CN performed most experiments on simulated and real datasets. AI designed the algorithm and performed the mathemathical proof of formula. SK implemented iLOCi program. AA designed the algorithm. CB performed the functional analysis of real dataset. PJS wrote the manuscript and discussed the test results. ST designed the algorithm, discussed the test results, and wrote the manuscript.

Electronic supplementary material

12864_2012_4412_MOESM1_ESM.pdf

Additional file 1: The mathematical details of ρdiff value and its relation with LD (iLOCi_details.pdf). This file includes the mathematical details of iLOCi formula and its relationship with the allele-based LD calculation. (PDF 372 KB)

12864_2012_4412_MOESM2_ESM.pdf

Additional file 2: Penetrance tables for dataset simulation (Penetrance_tables.pdf). This file includes the penetrance models used for dataset simulation of two-locus and high-order ineractions. (PDF 85 KB)

12864_2012_4412_MOESM3_ESM.xls

Additional file 3: Top 1000 SNP pairs from analyses of complete SNP set of WTCCC (TopPairs_Complete.xls). This file includes the list of top 1000 SNP pairs with their associated genes obtained from the iLOCi analyses of all SNPs passing the quality control step. The evidences for disease association of each identified gene as reported in WTCCC original paper and HuGE Navigator database are also shown. The genes identified as candidate disease genes from ToppGene prioritization are indicated with their rank numbers and P-values. (XLS 1 MB)

12864_2012_4412_MOESM4_ESM.xls

Additional file 4: Top 1000 SNP pairs from analyses of gene-only SNP set of WTCCC (TopPairs_GeneOnly.xls). This file includes the list of top 1000 SNP pairs with their associated genes obtained from the iLOCi analyses of gene-only SNPs. The evidences for disease association of each identified gene as reported in WTCCC original paper and HuGE Navigator database are also shown. The genes identified as candidate disease genes from ToppGene prioritization are indicated with their rank numbers and P-values. (XLS 2 MB)

12864_2012_4412_MOESM5_ESM.xls

Additional file 5: Pathway enrichment analysis of WTCCC datasets (Pathway_analysis.xls). This file includes the list of enriched biological pathways obtained from ToppGene program using the training sets of HuGE Navigator disease-associated genes. The pathway P-value is reported along with the list of iLOCi identified genes associated with such pathway. For each pathway, the number of genes previously reported in HuGE Navigator database, reported in WTCCC paper, and the novel disease genes, is shown. (XLS 175 KB)

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Piriyapongsa, J., Ngamphiw, C., Intarapanich, A. et al. iLOCi: a SNP interaction prioritization technique for detecting epistasis in genome-wide association studies. BMC Genomics 13 (Suppl 7), S2 (2012). https://doi.org/10.1186/1471-2164-13-S7-S2

Download citation

  • Published:

  • DOI: https://doi.org/10.1186/1471-2164-13-S7-S2

Keywords