Identifying gene-gene interaction is a hot topic in genome wide association studies. Two fundamental challenges are: (1) how to smartly identify combinations of variants that may be associated with the trait from astronomical number of all possible combinations; and (2) how to test epistatic interaction when all potential combinations are available. We developed AprioriGWAS, which brings two innovations. (1) Based on Apriori, a successful method in field of Frequent Itemset Mining (FIM) in which a pattern growth strategy is leveraged to effectively and accurately reduce search space, AprioriGWAS can efficiently identify genetically associated genotype patterns. (2) To test the hypotheses of epistasis, we adopt a new conditional permutation procedure to obtain reliable statistical inference of Pearson' s chi-square test for the contingency table generated by associated variants. By applying AprioriGWAS to age-related macular degeneration (AMD) data, we found that: (1) angiopoietin 1 (ANGPT1) and four retinal genes interact with Complement Factor H (CFH). (2) GO term “glycosaminoglycan biosynthetic process” was enriched in AMD interacting genes. The epistatic interactions newly found by AprioriGWAS on AMD data are likely true interactions, since genes interacting with CFH are retinal genes, and GO term enrichment also verified that interaction between glycosaminoglycans (GAGs) and CFH plays an important role in disease pathology of AMD. By applying AprioriGWAS on Bipolar disorder in WTCCC data, we found variants without marginal effect show significant interactions. For example, multiple-SNP genotype patterns inside gene GABRB2 and GRIA1 (AMPA subunit 1 receptor gene). AMPARs are found in many parts of the brain and are the most commonly found receptor in the nervous system. The GABRB2 mediates the fastest inhibitory synaptic transmission in the central nervous system. GRIA1 and GABRB2 are relevant to mental disorders supported by multiple evidences.
Gene-gene interactions have been proposed as one potential explanation of the well-known problem of missing heritability [1], and a recent report [2] has quantitatively demonstrated that possibility. Researchers have long attempted to identify interactions, with methods ranging from evolutionary genetic studies [3], [4], systems biology studies of model microbes [5] and quantitative genetic studies of inbred model organisms, to linkage [6] and association studies in human populations [7]–[14]. Although the definitions of the term “epistasis” used by biologists (Batson 1909) [15] and statisticians (Fisher 1918) [16] are different, they have the same consequences regarding different distributions of genotype patterns among different phenotypes. The main obstacle of interaction analysis is that the large number of multi-locus genotype combinations generated from large numbers of genetic variants is too high for current computational resources. This is in fact a well-known computational problem, known in the field of computer science as the ‘curse of dimensionality’ [17]. In this work we developed AprioriGWAS, a tool to address this problem. This tool is based on a successful algorithm in the field of computer science, Apriori [18]. Apriori was originally designed for supermarket data mining to assist shop owners in designing the layout of displayed products. Given customers' transactions, the algorithm can identify sets of items that frequently co-exist in transactions. For example, by knowing that customers usually buy milk and bread together, the shop owner can put them near each other in the store. Before describing the algorithm, we briefly give definitions of a few key terms: item is defined as an individual product, for example, bought by a customer; itemset stands for a set of items purchased together; length of itemset is defined as the number of items in the itemset. The process of growing a short itemset to a longer itemset is referred to as pattern growth. Generally, the key insights of Apriori are that: (1) frequent itemset with many items can be gained by growing itemset of short length; and (2) since subsets of any frequent itemset should also be frequent during pattern growth, itemsets predicted not to have any effect can be dropped during pattern growth, thereby significantly reducing the search space. In the case of GWAS, the number of individual genotypes is analogous to the number of transactions in supermarket data. The genotype of a variant is an item, and genotype combinations of different variants are an itemset, here also called a genotype pattern. Instead of just finding frequent genotype patterns, we want to find genotype patterns with different frequencies in cases and controls. We call them differential genotype patterns. While Apriori originally works on one database to find the most frequent itemsets, we are interested in patterns with different frequencies in two databases (cases and controls). To assess whether a pattern should be retained during pattern growth, we make use of the proportion test [19] (Methods). Interaction among variants is carried out after obtaining all differential genotype patterns. We test the possibility of interaction among variants involved in a differential genotype pattern by conducting Pearson' s Chi-square test for the contingency tables composed of all genotype patterns found for variants and phenotypes (Methods). In this step, we try to distinguish whether a differential pattern is caused by variants with marginal effects or by interaction effect. The process of pattern growth helps to narrow down the number of variant combinations to be tested for interaction effect. Using simulations following Marchini et al' s procedure [11], we demonstrate that AprioriGWAS can approximately achieve the same coverage of associated patterns as an exhaustive search, but with far lower CPU time. Determining all potential combinations that are statistically associated with disease does not automatically identify genuinely interacting genes. The daunting number of all combinations of variants heavily increases the load of multiple tests and mixes genuine signals with noise. As summarized by Anderson [20], in the regression model with two main effects terms and one interaction term, there is no exact permutation method for testing the significance of the interaction term. Buzkova et al [21] proposed a parametric bootstrap test for gene-gene and gene-environment interactions, which unfortunately is not practical for very large numbers of possible combinations of variants. Computer simulation [22] shows that whenever a trait is controlled by more than a single factor, it becomes possible for a neutral variant together with a major-effect variant as a pattern to be more strongly associated with the trait than with any of the causative factors [13]. These indirect associations are true associations for statistical purposes, and can be indistinguishable from medical causative associations [22]. To distinguish general association and interaction effects, we developed a new conditional permutation test to distinguish genuine interactions from the artifacts generated by the combination of a major-effect variant with a neutral variant (Methods). We demonstrate that our new approach has a magnitude lower false discovery rate (FDR) compared with regular permutation, while maintaining comparable power. We applied AprioriGWAS to age-related macular degeneration (AMD [MIM 153800]), which has been deemed a good example of a small number of common variants explaining a large proportion of heritability [1]. Among the most significant patterns, we found six pairs of retinal genes interacting with each other. An exciting example is the interaction of a gene involved in an AMD treatment target, ANGPT1, with another important AMD gene, CFH. Overall, the potentially interacting genes were enriched in glycosaminoglycan biosynthetic process (). Many studies have shown that the interaction between glycosaminoglycans (GAGs) and CFH plays an important role in the disease pathology of AMD. We also applied AprioriGWAS to bipolar disorder; we found potential interactions inside individual gene (8 out of 18 genes are related with mental disorder) and interactions across gene or chromosomes. Further results will be presented in full later. The remainder of this paper is organized as follows. In the next section we introduce the AprioriGWAS algorithm for mining possible interaction variants, as well as the conditional permutation approach for testing interactions. We then evaluate the performance of AprioriGWAS with simulated data and compare it with logistic regression implement in Epistasis function of PLINK. Lastly we demonstrate applications of AprioriGWAS to AMD and WTCCC bipolar data and exciting findings from both datasets.
Historically, the Apriori algorithm can be traced back to the seminal paper published by IBM Research in 1993 [18]. The concept of the main technique is that a subset of frequent itemset should also be frequent. Based on this concept, frequent itemset with more items may be found by stepwise growth of smaller frequent itemset, which saves substantial computational resources. Interested readers may refer to their original paper [18] for a professional description or to our own longer report [23] for illustrative descriptions. Here we briefly outline the main steps. Suppose one wants to mine frequent itemset with length no more than n. Apriori will usually scan dataset in n rounds (unless there is no new frequent itemset generated in a certain round before n, thereby forcing the algorithm to halt). In the first round, it will initiate the 1-itemsets that are frequent. In each subsequent round, it will take the frequent itemset generated in the last round as starting point and grow any itemset by adding one more item. Retention of the new itemset will be decided by firstly predicting how likely it will be and then, given a positive prediction, by checking the actually supporting transactions. Finally, the collection of all frequent itemset in all rounds will be reported. In this paper, genotype patterns are defined as genotype combinations of different variants. We use integer numbers as ids of variants; then we can have, for instance, a pattern like 46AT_609GG_1099CC, denoting a pattern composed of a variant with id 46 and genotype AT combined, a variant with id 609 and genotype GG, and a variant with id 1099 and genotype CC. The key goal is to find genotype patterns that have a significant frequency difference in cases and controls (called differential patterns in this paper). The algorithm of AprioriGWAS is divided into two steps. First, detecting differential genotype patterns by an Apriori-like strategy. Obviously, the same set of variants can lead to several differential genotype patterns. Second, testing interaction among a set of variants by testing association of all possible combinations of genotype patterns against case/control status. The first step helps to narrow down the combinations of variants need to be tested. Due to multiple test problems and potential association of single variants involved in the differential genotype pattern, we adopt a new conditional permutation in the second step to control the marginal effect of single variants for testing of variant interactions.
We simulated data by two-locus interaction models proposed by Marchini et al [11] (Methods), in which three types of interactions are generated. We then applied regular permutation and conditional permutation to control family-wise type I error. The performances of regular permutation and conditional permutation test (Methods) are demonstrated in Figure 1A and 1B. We compared both power and FDR, using regular permutation and conditional permutations to adjust thresholds for type I error. Family-wise type I error was set to 0. 05 for both methods. It is evident that the FDR was significantly reduced by the conditional permutation test, although some power is sacrificed compared with regular permutation. To demonstrate that the nominal p-value of a contingency table for multi-variants could be in large part caused by individual variants with strong marginal effect, we took a real example from analyzed AMD data. Figure 2A shows two variants, each with no marginal effect, but in combination with strong marginal effect. Figure 2B shows two variants, one has strong marginal effect, and the other does not show any marginal effect. Although the nominal p-value of the contingency table is more significant than the pair of variants in Figure 2A, one can deduce that the low p-value from Figure 2B is in large part caused by the variants with strong marginal effect; in Figure 2A, on the other hand, there must be some interaction effect. As mentioned, AprioriGWAS manages to dramatically speed up the search process by dropping the candidate genotype patterns unlikely to grow to differential pattern. Since it is based on prediction at an early stage in the search, it still theoretically runs the risk of mistakenly dropping sensible patterns. Here we quantitatively tested the percentage of mistakenly dropped differential patterns by comparing AprioriGWAS and exhaustive search (Method). Figure 3 shows the comparison between searching for combinations of variants (with default parameters in AprioriGWAS) and exhaustive search. We found that 97% of all differential genotype patterns found by exhaustive search were covered by the results from AprioriGWAS. With such high coverage, the chance of losing possible interaction variants is minimized. There are a few points below 85%, reflecting that there is variation of power to cover all potential combinations. It is true that the overall coverage is subject to lots of parameters, like sample size and allele frequency. To minimize this variation, larger sample size is always desirable. We compared the ability of AprioriGWAS to find interacting variants with traditional single locus genotypic test and exhaustive search in PLINK [28] (epistasis function). The epistasis function in PLINK for case control data is basically stepwise logistic regression. We chose to use the all combinations option. The power comparison is based on two levels: finding at least one casual variant, or finding both interacting variants (Figure 4). For Level 1, detecting at least one causal variant, we found that the traditional single variant test had the highest power in Model 1, which has explicit marginal effects for both causal variants. AprioriGWAS performed similarly with the single loci test in Model 2, and had better power in Model 3 (Figure 5). This is natural, since Model 2 and 3, which contain no explicit marginal effects, are expected to be harder to detect without an interaction-based searching strategy. For Level 2, detecting both interacting variants, it is evident that AprioriGWAS had the highest power in most cases of Model 2 and 3 (Figure 4). On the other hand, the performance of the epistasis function in PLINK, which exhaustively searches all combinations, was not as good in all cases. This is because: (1) stepwise logistic regression does not capture the interactions well, since the effects of the terms are added in a linear manner, whereas AprioriGWAS explicitly addresses detailed patterns; (2) in stepwise logistic regression the genuine interactions are buried by the noise of a too large number of combinations, whereas with the conditional permutation test used in AprioriGWAS, genuine interactions are able to stand out. When comparing corresponding panels in Figure 4 and Figure 5, it is observed that for the single variant test the power of finding both interacting variants (i. e. , Level 2) dropped significantly compared with the power of finding at least one causal variant (i. e. , Level 1). By contrast, interaction based methods, i. e. , both AprioriGWAS and PLINK epitasis, maintained similar power for both levels. This was not unexpected since the interaction-based strategies should be better able to find an epistasis effect. We also simulated data that have more SNPs (1,000,000) and find that the relative power between three methods and interaction models remain similar although the absolute powers are all decreased. (Figure S1) Figure 6 shows the power of AprioriGWAS and single variant test on three classical genetic models studied in model organisms. There are three powers for each genetic models: power for detecting at least one gene using single variant test, power for detecting both genes using single variant test, and power for detecting both genes using AprioriGWAS. Since PLINK is not scalable for such a dataset, we have not achieved power estimates for logistic regression. For the model “Duplicated Dominant”, AprioriGWAS outperforms single marker test for detecting single gene or both genes, whereas for models “Duplicated Recessive” and “Dominant & Recessive Interaction”, AprioriGWAS is more powerful for detecting both genes, but not for detecting single genes. It is notable that the power of detecting both genes in the model “Dominant & Recessive Interaction”, in which epistasis is functioning; single variant test has almost zero power (0. 1%) while AprioriGWAS has around 50% power. We compared the speed of our method with the epistasis function in PLINK. Figure 7 shows that the default threshold setting in AprioriGWAS was approximately a magnitude faster. Although retaining candidate genotype patterns in memory can help speed up the algorithm, its affordability is subject to the particular computational resources. We took the strategy of writing candidate patterns on hard disk for each round of pattern extension. The genotype data used to be relatively small comparing with the patterns however is getting larger and larger empowered by new sequencing platforms. To solve this problem, we implemented AprioriGWAS using HDF5-based data format [29] which stores genotype data on disk and accesses them as though stored in main memory. Therefore, the memory usage is scalable to whatever size of potential dataset and the speed is not scarified. (See more on computational and memory complexity in section Discussion.) We applied AprioriGWAS on published AMD data [26]. We identified 168 significant pairs of variants (family-wise type I error = 0. 01), presented in Table S1. By checking published functional literals and gene annotations, as well as GO enrichment of the genotype patterns, we learned that the findings are well validated by existing functional studies and clinical applications. Besides AMD data that were extensively analyzed by the community interested in gene-gene interactions, we also applied AprioriGWAS on Bipolar Disorder data from WTCCC [27] to further test whether it is scalable for larger dataset. The whole task was distributed onto 1,000 CPUs in a cluster and the average execution time for a single job is 56. 8 hours. Only 4 Gb memories were employed during the computation, evidencing the great performance of HDF5-based implementations.
Regardless the goal being interaction or single gene, statistical tests all suffer from the problem of false positives. Since the numbers of variants (and their combinations) are usually a few magnitudes larger than the sample size for most association studies, it will be common to see false positives. The current practice in the community is that researchers who would like to claim association or carry out experimental validations usually have to check whether the results are replicable in other independent dataset (s) Researchers who use AprioriGWAS can also use this to filter results before doing experimental validations. As an example, we use another independent dataset for AMD study [64] to check whether the results are replicable. Among the five interactions with CFH reported in this paper, we found that BBS9/CFH and CHRM2/CFH are replicated in the other dataset. However, we understand that these two datasets are very different: one is wet AMD and the other is dry AMD. One of them is more prevalent in Asia than the other. Therefore, our further analysis of data in [64] may not serve as perfect replication of the findings presented, although it suggests that BBS9 and CHRM2 may be of higher priority for further experimental validations. The most commonly used multiple variants analysis is stepwise regression, in which variants are added to the regression equation one after another by some suitable criteria. But statistical analysis shows that the usual stepwise model selection methods are path dependent and therefore suboptimal [65]. Besides regression, some methods are based on discrete mathematics, like the Combinatorial Partition Method (CPM) [66] and its refined version, the Restricted Partition Method (RPM) [67]. However, RPM still requires a daunting number of tests when the number of variants is high. This is because its insight into reducing tests lies in its practice to combine close phenotypes, which consequently does not entirely solve the problem of too many combinations of genotypes. Another well-known method of counting potential combinations is multifactor dimensionality reduction (MDR). It collapses cells in a contingency table into two groups and conducts a test on them. Essentially however it reduces the dimensionality of testing, rather than reducing the dimensionality of the process of counting genotype patterns. Therefore, when the number of variants is large, it still suffers from the “curse of dimensionality” [17]. Bayesian methods leveraging MCMC, e. g, BEAM [53] or epiMODE [8], should theoretically suffer less from computational limitations, but they do not directly test detailed combinations of genotype patterns and thereby sacrifice the advantages of fine scale learning of gene-gene interactions. Another branch of frequently used methods is two-stage analysis [68], by which the investigator can utilize relatively “simple” or computationally efficient tests to choose qualified variants in the first stage analysis. Then, taking advantage of the relatively small number of variants, the investigator can adopt some advanced but computationally heavy test to identify interacting genes. However, due to a lack of strong prior knowledge, the true signals might have been removed from the first stage if the procedure was not well designed. As an example, interacting variants with no marginal effect may be filtered out if one uses tests based on marginal effects of single variants in the first stage. Nevertheless, with good design, this approach is still very promising and can be combined with all the approaches reviewed above; and it can naturally also be combined with the method proposed in this work. Computation time and spatial complexities of the tool may be interesting to the reader. The number of transactions for original Apriori corresponds to sample size in GWAS; the number of items is equivalent to the number of variants and the itemsets. In contrast to supermarket data, GWAS data have a limited number of “transactions”, but a large number of “items” in two datasets, cases and controls. Both conditions make the problem more difficult. The time spent reading the data in each round of pattern growth is constant. In addition, the computational resources cost depends on how many combinations of genetic variants will be generated and tested. The more combinations are tested, the less likely it is that genuine patterns are missed, though of course more resources will be used. In AprioriGWAS, there are several parameters for the user to specify according to their computer resources and understanding of the disease model. The threshold for the proportion test and minimal support of concerned itemsets are parameters that affect candidate search space, algorithm speed, and power of detecting all distinct genotype patterns. When these parameters are set to zero, AprioriGWAS will exhaustively search all possible combinations. (Please refer to our Manual of AprioriGWAS for the tradeoffs and discussions on setting these parameters according to computational resources.) Those familiar with Apriori may suggest that, given Apriori' s ability to also mine association rules, one could also treat the case control label as items and directly adopt Apriori for case/control data. The result will then be a subset of variants that can imply the case/control labels. But searching frequent itemsets and then mining the association between genotype pattern and disease status is inefficient, since frequent genotype patterns are not necessarily associated with phenotype; on the other hand, genotype patterns strongly associated with phenotype may not necessarily be in high frequency, and such an association could be distributed in different patterns than the same variants combinations. Instead of the conditional permutation proposed here, one could also consider Bonferroni correction. For n variants with search length of m, the total number of combinations is huge. Given the natural correlation of the combinations, it is clearly far more stringent than necessary. However, only correcting on the number of differential pattern tested produces a bias in the other direction, since the nominal value of the significance level of the chi-square test for the contingency table will be inflated by the selection procedure [69]. It is therefore always preferable to use a permutation test for the whole procedure. With regular permutation, one permutes the Case/Control label and then performs the whole test process. The smallest P-value of each permutation are ranked, allowing one to get the distribution of test statistics under “Null” from the permuted dataset. With regular permutation, no variant should have marginal effect, and the p-value of the contingency table for the combination of variants is under the null hypothesis of no variants having marginal effect. However, regular permutation suffers from an inflated significance level for contingency tables containing variants with marginal effects. This is due to the fact that when a contingency table is composed of at least one variant with strong marginal effect, the p-value for that contingency table becomes extremely small compared with regular permutation results. The FDR is therefore very high, even close to 1. To solve the problem of an inflated significance level by a contingency table composed of at least one variant, v, with strong marginal effect, we developed a conditional permutation procedure (Methods), which helps get the null distribution of the p-value of a contingency table composed of the variant and other variants. Simulation results show that, when we control the family-wise type I error by conditional permutation, we also keep FDR well controlled. Compared with INTERSNP [13], which lists only the top 50 variant combinations including the variant with marginal effect, conditional permutation in AprioriGWAS keeps FDR well controlled in a systematic way. Another concern might be whether these differential genotype patterns are artifacts caused by linkage disequilibrium (LD). We believe this is not the case, since the LD should impact both cases and controls, and therefore the pattern created by LD will not be differential unless the LD structure is significantly different in cases and controls for particular genetic variants. If that is the case, then there must be some reason of selection to explain the deviation in the genotype pattern, and it is difficult to judge whether this is an artifact or something of interest. In addition, our conditional permutation also breaks LD between interacting variants. Low-frequency or rare variation might impact the performance of the method, even when explicitly only testing for interactions among common variants. What matters is the extent of LD between causal rare variants and testing common variants. We haven' t addressed this problem in the current method. It would be interesting to extend AprioriGWAS toward that direction. There may be non-trivial statistical challenges since the low-frequency or rare variants are usually less shared by the individuals therefore their combinations that form genotype patterns will be even less shared by individuals. For a given set of variants, we will have many patterns with little supports.
Genes do not operate in vacuum. They interact with each other in many ways. Therefore, to figure out genetic causes of disease by case-control association studies, it is important to take interactions into account. There are two fundamental challenges in interaction-focused analysis. The first is the number of possible combinations of genetic variants easily goes to astronomic which is beyond current computational facility, which is referred as “the curse of dimensionality” in field of computer science. The other is, even if all potential combinations could be exhaustively checked, genuine signals are likely to be buried by false positives that are composed of single variant with large main effect and some other irrelevant variant. In this work, we propose AprioriGWAS that employees Apriori, an algorithm that pioneers the branch of “Frequent Itemset Mining” in computer science to cope with daunting numbers of combinations, and conditional permutation, to enable real signals standing out. By applying AprioriGWAS to age-related macular degeneration (AMD) data and bipolar disorder (BD) in WTCCC data, we found interesting interactions between sensible genes in terms of disease. Consequently, AprioriGWAS could be a good tool to find epistasis interaction from GWA data.
algorithms
genome-wide association studies
computer and information sciences
mathematics
statistics (mathematics)
gene regulatory networks
genome analysis
genetics
applied mathematics
biology and life sciences
physical sciences
computational biology
biostatistics
genetics of disease
human genetics
AprioriGWAS, a New Pattern Mining Strategy for Detecting Genetic Variants Associated with Disease through Interaction Effects