Detecting purely epistatic multilocus interactions by an omnibus permutation test on ensembles of twolocus analyses
 Waranyu Wongseree^{1},
 Anunchai Assawamakin^{2},
 Theera Piroonratana^{1},
 Saravudh Sinsomros^{1},
 Chanin Limwongse^{2} and
 Nachol Chaiyaratana^{1, 2}Email author
DOI: 10.1186/1471210510294
© Wongseree et al; licensee BioMed Central Ltd. 2009
Received: 5 March 2009
Accepted: 17 September 2009
Published: 17 September 2009
Abstract
Background
Purely epistatic multilocus interactions cannot generally be detected via singlelocus analysis in casecontrol studies of complex diseases. Recently, many twolocus and multilocus analysis techniques have been shown to be promising for the epistasis detection. However, exhaustive multilocus analysis requires prohibitively large computational efforts when problems involve largescale or genomewide data. Furthermore, there is no explicit proof that a combination of multiple twolocus analyses can lead to the correct identification of multilocus interactions.
Results
The proposed 2LOmb algorithm performs an omnibus permutation test on ensembles of twolocus analyses. The algorithm consists of four main steps: twolocus analysis, a permutation test, global pvalue determination and a progressive search for the best ensemble. 2LOmb is benchmarked against an exhaustive twolocus analysis technique, a set association approach, a correlationbased feature selection (CFS) technique and a tuned ReliefF (TuRF) technique. The simulation results indicate that 2LOmb produces a low falsepositive error. Moreover, 2LOmb has the best performance in terms of an ability to identify all causative single nucleotide polymorphisms (SNPs) and a low number of output SNPs in purely epistatic two, three and fourlocus interaction problems. The interaction models constructed from the 2LOmb outputs via a multifactor dimensionality reduction (MDR) method are also included for the confirmation of epistasis detection. 2LOmb is subsequently applied to a type 2 diabetes mellitus (T2D) data set, which is obtained as a part of the UK genomewide genetic epidemiology study by the Wellcome Trust Case Control Consortium (WTCCC). After primarily screening for SNPs that locate within or near 372 candidate genes and exhibit no marginal singlelocus effects, the T2D data set is reduced to 7,065 SNPs from 370 genes. The 2LOmb search in the reduced T2D data reveals that four intronic SNPs in PGM1 (phosphoglucomutase 1), two intronic SNPs in LMX1A (LIM homeobox transcription factor 1, alpha), two intronic SNPs in PARK2 (Parkinson disease (autosomal recessive, juvenile) 2, parkin) and three intronic SNPs in GYS2 (glycogen synthase 2 (liver)) are associated with the disease. The 2LOmb result suggests that there is no interaction between each pair of the identified genes that can be described by purely epistatic twolocus interaction models. Moreover, there are no interactions between these four genes that can be described by purely epistatic multilocus interaction models with marginal twolocus effects. The findings provide an alternative explanation for the aetiology of T2D in a UK population.
Conclusion
An omnibus permutation test on ensembles of twolocus analyses can detect purely epistatic multilocus interactions with marginal twolocus effects. The study also reveals that SNPs from largescale or genomewide casecontrol data which are discarded after singlelocus analysis detects no association can still be useful for genetic epidemiology studies.
Background
Complex diseases cannot generally be explained by Mendelian inheritance [1] because they are influenced by genegene and geneenvironment interactions. Many common diseases such as asthma, cancer, diabetes, hypertension and obesity are widely accepted and acknowledged to be results of complex interactions between multiple genetic factors [2]. Attempts to identify factors that could be the causes of complex diseases have led to many genomewide association studies [3, 4]. Raw results from these attempts produce a large amount of single nucleotide polymorphism (SNP) data from every individual participating in the trials.
For genetic epidemiologists, data sets from genomewide association studies present many challenges, particularly the correct identification of SNPs that associate with the disease of interest from all available SNPs [5]. This challenge can be treated as a pattern recognition problem which aims to identify an attribute or SNP set that can lead to the correct classification of recruited samples. Heidema et al. [5] and Motsinger et al. [6] have reviewed and identified many machine learning techniques that are suitable to the task. Among many strategies and techniques, the protocol that appears to be most promising for genomewide association studies involves two main steps: SNP set reduction and classification model construction [7]. From a machine learning viewpoint, attribute selection techniques can be divided into three main categories: filter, wrapper and embedded approaches [8]. In a filter approach, a measure or an index is used to determine the correlation between attributes and classes, e.g. affected and unaffected status in a casecontrol study. Attributes that are deemed to be important for the classification according to the measure are then selected. The filter approach includes χ^{2} and odds ratio tests [9, 10], omnibus permutation tests [11–13], a correlationbased feature selection technique [14], a ReliefF technique [15] and a tuned ReliefF technique [16]. In a wrapper approach, the significance of an attribute subset is evaluated from the classification performance by a classifier. The capability of the wrapper approach to identify significant attributes thus depends on the chosen classifier and the search algorithm for the identification of the best attribute subset. Combinatorial [17] and restricted partitioning methods [18], a multifactor dimensionality reduction method [19–25] and a polymorphism interaction analysis technique [26] are examples of the wrapper approach. An embedded approach concentrates on informative attributes during the construction of a classification model. Examples of the embedded approach include a genetic algorithm with Boolean algebra [27], genetic programming based decision trees [28, 29], random forests [30–32] and evolutionary neural networks [33, 34]. Based on this categorisation, classification models are not direct outputs from filterbased techniques. On the other hand, classification models are readily prepared as outputs from the wrapper and embedded approaches. In other words, the last two approaches can also be regarded as classification model construction techniques.
The success of the twostep pattern recognition approach relies heavily on the attribute selection step [14]. In casecontrol studies, epistatic effects play a vital role in establishing the difficulty level of SNP screening problems [35, 36]. Epistasis in the simplest form can be represented by disease models that require genotype inputs from two interacting SNPs [37, 38]. Many attempts have been made to produce consistent definitions and categorisation of different types of epistasis models [2, 35, 39–41]. According to Musani et al. [2], a pure epistasis model [42] is difficult because each SNP exhibits no marginal singlelocus effect in the model. As a result, it is impossible to detect the pure epistasis by univariate statistical tests. Examples of complex diseases that casecontrol studies have uncovered putatively pure epistasis include type 2 diabetes mellitus (T2D) [43–46] and metabolic syndrome [47]. Due to the difficulty of screening for each SNP independently, it is suggested that attention should be focused on the analysis of differences between twolocus genotype distribution within case and control groups [40] and multilocus Bayesian statistical analysis [48, 49].
A number of SNP screening and association detection techniques have adopted the twolocus genotype monitoring strategy as their core engines [40, 50–52]. The search for interactions can be carried out via either exhaustive analysis [52] or the analysis that can be divided into two stages, incorporating singlelocus analysis for the prescreening purpose [40, 50, 51]. In the twostage mode, at least one SNP that involves in the construction of twolocus genotype unit must be a strong candidate for the association explanation, usually verified through univariate statistical tests. Each mode of the twolocus analysis possesses different strengths and weaknesses. The exhaustive analysis has a full capability of detecting pure epistasis but requires larger computational efforts [52]. In contrast, the twostage analysis is more practical for largescale data but with some risk of missing possible pure epistasis [50]. More practical usage of both twolocus analysis modes in real casecontrol studies is required before the feasibility issue can be fully addressed.
Many genetic association studies reveal that various complex diseases are results of putative multilocus interactions [11, 46, 53]. With the constraints on a computational capability, exhaustive multilocus analysis in largescale or genomewide association studies would be infeasible [52]. On the other hand, singlelocus analysis would be unsuitable for the detection of pure epistasis. One possible approach that provides a tradeoff between a computational limitation and an epistasis detection capability is to capture a multilocus interaction by combining multiple results from twolocus analysis. To achieve this, it is necessary to prove that once a multilocus interaction model is broken down into a combination of twolocus models, all or some of these models remain detectable through twolocus analysis. Although it is hinted in an early work on twolocus analysis [52] that the proposed approach is plausible, explicit experimentation and testing has never been conducted.
In this article, the feasibility of employing an ensemble of twolocus analyses for the multilocus interaction determination is demonstrated. Specifically, the significance of the twolocus analysis ensemble is assessed by an omnibus permutation test [54]. The proposed method is inspired by a set association approach [11], in which a limited number of sets that contain different numbers of SNPs are explored for possible association. These SNP sets are crucial in the global pvalue calculation of the selected set via a permutation test and thus the decision to accept or reject the null hypothesis of no association. In other words, SNP set exploration and selection is required to assess the significance of the identified association. This means that the set association approach is equally interested in both SNP set selection and testing for significant association. The primary function of the proposed method is to detect possible association and assess its significance through the exploration of different ensembles of twolocus analyses. Hence, the proposed method is also equally interested in both ensemble selection and testing for significant association.
The proposed method is benchmarked against a simple exhaustive twolocus analysis technique, the set association approach [11], the correlationbased feature selection technique [14] and the tuned ReliefF technique [16]. These filterbased attribute selection techniques are suitable for the benchmark trial since they are capable of detecting association. The casecontrol classification models constructed from screened SNPs via a multifactor dimensionality reduction method [19] are also provided.
Results and discussion
Algorithm
Twolocus analysis
Consider a casecontrol genetic association study with n_{ m }SNPs, for each pair of SNPs, a 2 × 9 contingency table with rows for disease status and columns for genotype configurations is created. A χ^{2} test statistic and the corresponding pvalue can subsequently be computed. With the total of n_{ m }SNPs, there are = n_{ m }!/((n_{ m } 2)!2!) possible SNP pairs. As a result, the pvalue from each twolocus analysis must be adjusted by a Bonferroni correction. The Bonferronicorrected pvalue from each analysis is the lower value between × the uncorrected pvalue and one.
Permutation test
where t is the number of permutation replicates which is set to 10,000 in this study and · denotes the size of a set.
Global pvalue determination
Search for the best ensemble of twolocus analyses
A simple progressive search is used to identify the best ensemble of twolocus analyses. The search begins by locating the best twoSNP unit with the smallest Bonferronicorrected χ^{2}'s pvalue from step 1. A permutation test is then performed for this twolocus analysis, yielding both raw and global pvalues since only one hypothesis has been explored. Next, the search attempts to combine the existing best twoSNP unit with the twoSNP unit that possesses the next smallest Bonferronicorrected χ^{2}'s pvalue from step 1 and does not have a higher permutation pvalue than the first twoSNP unit. If this new ensemble yields either a higher raw pvalue or the same raw pvalue but a higher global pvalue from a permutation test, the search is terminated and the association is explained by the previously identified twolocus analysis. Otherwise, the best ensemble of twolocus analyses is updated and the process of appending more twoSNP units to the ensemble continues. The progressive search terminates when deterioration in the raw or global pvalue is detected, or all possible twolocus analyses have been included in the ensemble. It is recalled from step 3 that for the best ensemble containing E  1 < twolocus analyses, its global pvalue is obtained from the evaluation of E hypotheses.
Validity of the algorithm
A permutation replicate in 2LOmb is constructed by randomly assigning the case or control status to each sample while maintaining the original proportion of case and control samples. Once the construction of a permutation replicate is finished, the assigned case and control labels remain fixed to the samples. The pattern of case and control labels in each permutation replicate is thus constant and unique. Therefore, the Bonferronicorrected χ^{2}'s pvalues from any twoSNP units within a permutation replicate are calculated from the same casecontrol data set. Hence, the combining of these Bonferronicorrected χ^{2}'s pvalues via a Fisher's combining function is attainable. The calculation of Fisher's test statistics from all permutation replicates and the original data set leads to the raw or unadjusted pvalue for the null hypothesis of the ensemble e as given in equation 1. Since the same set of permutation replicates is always used during the evaluation of each ensemble, the raw pvalues for the null hypotheses from all ensembles can be directly compared against one another. Furthermore, the global pvalue calculation is based on this set of permutation replicates. This is possible because the unadjusted pvalue for the permutation replicate i of ensemble e or can be calculated in a similar manner to the raw pvalue as defined in equation 2. The unadjusted pvalues for the same permutation replicate but different ensembles can also be directly compared and the subsequent calculation of is attainable. With and , the pvalue for the global null hypothesis that incorporates all E explored hypotheses can be determined by equation 3. In summary, only one set of permutation replicates is required for the calculation of both the raw pvalue for the null hypothesis of every ensemble and the global pvalue. The pvalues can be compared in each step of 2LOmb. Consequently, the selection of the best ensemble for association explanation can be carried out via a pvalue comparison.
Testing with simulated data
2LOmb is benchmarked against a simple exhaustive twolocus analysis technique, a set association approach (SAA) [11], a correlationbased feature selection (CFS) technique [14] and a tuned ReliefF (TuRF) technique [16] in a simulation trial. The exhaustive twolocus analysis is simply the twolocus analysis procedure from the first step of the 2LOmb algorithm. An interaction is declared if at least one twoSNP unit with a Bonferronicorrected χ^{2}'s pvalue below 0.05 is detected. The exhaustive twolocus analysis reports all SNPs that meet this detection condition. The simulation covers two main data categories: null data of no significant genetic association and data with causative SNPs which signify pure epistasis. The algorithm performance on the null data provides an indication for the falsepositive error. On the other hand, the algorithm performance on the data with causative SNPs indicates the detection capability. An efficient algorithm should produce an output with a low number of SNPs and a high number of correctlyidentified causative SNPs when epistasis is present. Similarly, it should also report that there are no causative SNPs in the null data. These two measures on the number of SNPs in the results are used as the performance indicators.
Each simulated data set contains 1,000, 2,000 or 4,000 SNPs in which either there are no causative SNPs or there is pure epistasis, governed by two, three or four causative SNPs. The allele frequencies of all causative SNPs are 0.5 while the minor allele frequencies of the remaining SNPs are between 0.05 and 0.5. The data set consists of balanced casecontrol samples of sizes 400, 800 or 1,600. All SNPs in control samples are in HardyWeinberg equilibrium (HWE) [59]. The genotype distribution of causative interacting SNPs follows the pure epistasis model by Culverhouse et al. [42], leading to three interesting values of heritability: 0.01, 0.025 and 0.05. Every SNP in each data set exhibits no marginal singlelocus effect (Bonferronicorrected χ^{2}'s pvalue > 0.05). Twentyfive independent data sets for each simulation setting are generated via a genomeSIM package [60]. A paired ttest is suitable to assess the significance of results since the same simulated data sets are used during the algorithm benchmarking.
The performance of many existing attribute selection techniques for pattern recognition depends on the level of attribute interactions. A number of techniques, including CFS, appear to function well under a moderate level of interactions. However, the performance of CFS appears to be significantly reduced when the interaction level becomes too high [14, 61] because CFS favours an attribute that is strongly correlated with the classification outcomedisease status in this studywhile at the same time is not correlated with other attributes. Since the main driving force behind epistasis is the interaction between SNPs, which are themselves attributes, CFS would not intuitively select all causative SNPs. Consequently, the SNP set produced by CFS appears to contain only uncorrelated SNPs. Obviously, a SNP that is a part of the interaction model would occasionally be picked up by CFS but CFS never successfully identifies all causative SNPs in any interaction problems. In addition, CFS reports more erroneous SNPs than other techniques in the null data problem and all three interaction problems due to many SNPs being uncorrelated.
The benchmarking of attribute selection techniques by Hall and Holmes [14] also reveals that ReliefF [15] is better than CFS in problems with a high level of interactions. Since ReliefF is essentially the core engine of TuRF, the results from this study are in agreement with the early benchmark trial. This finding strengthens the observation that the interaction level of SNPs in pure epistasis models is too high for CFS to handle. Similar to its predecessor, the performance of TuRF still depends on both the number of attributes and sample size. TuRF performs well in the majority of simulation scenarios with 1,0002,000 SNPs and 8001,600 samples. These scenarios are relatively easy since the number of SNPs is small while the sample size is large. However, the size of output SNP set, reported by TuRF from the null data problem and all three interaction problems, increases significantly when the difficulty level rises by either reducing the sample size or increasing the number of SNPs. This implies that when the problem contains a large number of candidate SNPs, the only way to ensure that TuRF reports a proper SNP set is to use a relatively large sample size, making it impractical in real genetic association studies due to many factors including disease prevalence, population size and genotyping cost.
The global pvalues in most of the SAA results from the null data problem and all three interaction problems exceed 0.05, showing that SAA reports a low falsepositive result in the null data problem. Nonetheless, SAA remains unsuitable for detecting pure epistasis because of its high falsenegative error. This poor performance can be traced back to the manner in which SAA exploits an omnibus permutation test. As stated earlier, singlelocus analysis does not detect any association between a SNP and the disease in this study. Hoh et al. [11] have demonstrated that genetic association can be more significantly observed when the singlelocus test statistics are combined together. Nonetheless, there is an additional requirement that each causative SNP must exhibit a marginal singlelocus effect. In the current study, the association signal from each causative SNP is lower than the required threshold, leading to similar test statistics and global pvalues for both combinations of multiple SNPs which include causative SNPs and those which exclude causative SNPs.
Both 2LOmb and exhaustive twolocus analysis technique are capable of identifying all causative SNPs. However, the size of output SNP set from 2LOmb is significantly smaller than that from the exhaustive twolocus analysis. Appended SNPs to the causative SNPs in the output from 2LOmb and those from the exhaustive twolocus analysis are erroneous SNPs. These erroneous SNPs are parts of false twoSNP units with Bonferronicorrected χ^{2}'s pvalues less than 0.05. A similar trend of results regarding the size of output SNP set is also observed in the benchmark trial involving the application of 2LOmb and exhaustive twolocus analysis to the null data. This signifies that the permutation test and the progressive search embedded in 2LOmb can help reducing the number of erroneous SNPs in the output.
As mentioned earlier, 2LOmb produces the best results among five techniques in the benchmark trial. 2LOmb has a low falsepositive error in the null data problem and is capable of detecting all causative SNPs in every simulated data set in all three interaction problems. This performance is further strengthened by highly significant global pvalues in 2LOmb results from all three interaction problems (p < 0.0001) and the presence of a SNP in common among some or all pairs of twoSNP units in the three and fourlocus interaction problems. Nonetheless, some of the 2LOmb outputs contain a few erroneous SNPs which are irrelevant to the correct association explanation. Since all three interaction problems involving different numbers of causative SNPs are investigated by varying the total number of SNPs, the sample size and the level of heritability, these parameters may influence the number of erroneous SNPs in the 2LOmb results. Similarly, the total number of SNPs and the sample size may affect the number of erroneous SNPs in the 2LOmb results from the null data problem. ANOVA reveals that the only source of variation that significantly affects the number of erroneous SNPs in the null data, twolocus interaction and threelocus interaction problems is the sample size (p < 0.000001). In addition, the sample size must be greater than 800 for an increase in the number of erroneous SNPs to be significant. In contrast, ANOVA reveals that two sources of variation that affect the number of erroneous SNPs in the fourlocus interaction problem are the sample size (p < 0.000001) and the total number of SNPs (p < 0.00005). Similar to the null data, twolocus interaction and threelocus interaction problems, the sample size in the fourlocus interaction problem must be greater than 800 to create a significant increase in the number of erroneous SNPs. On the other hand, the number of erroneous SNPs appears to decrease when the total number of SNPs increases. These two sources of variation also interact with each other (p < 0.005). However, the interaction is most evident only when the sample size is large, i.e. when the sample size is 1,600.
ANOVA shows that the number of erroneous SNPs in the 2LOmb results is influenced by the sample size and the total number of SNPs but not by the heritability. It is observed that the number of erroneous SNPs increases when the sample size is large. This counterintuitive phenomenon can be explained as follows. As 2LOmb combines pvalues that are determined from χ^{2} tests, the number of entries for the contingency table construction is large when the sample size is large. This subsequently leads to a significantly large χ^{2} statistic and hence an extremely small pvalue if the SNPs under consideration are causative SNPs. At the same time, the possibility that a reasonably large χ^{2} statistic and a small pvalue can be obtained by chance from a twoSNP unit which is irrelevant to the correct association explanation also inevitably increases. With the increase in the possibility of erroneous SNP inclusion, the size of output SNP set gets bigger when the sample size is large. Another observation that appears to be counterintuitive is the reduction in the number of erroneous SNPs when the total number of SNPs increases. This phenomenon is the result of the Bonferroni correction usage. When the total number of SNPs is doubled, the Bonferroni correction factor in 2LOmb is quadrupled. A higher correction factor leads to a more stringent criterion for SNP selection. This subsequently leads to the reduction in the number of erroneous SNPs when the total number of SNPs is large.
In contrast to the first two parameters, different levels of heritability appear to have no effect on the 2LOmb results because all simulated data sets have balanced casecontrol samples and the embedded interaction models have the same architecture. For instance, a twolocus interaction model leads to zero penetrances for genotypes AABB, AABb, AaBB, Aabb, aaBb and aabb. Hence, the penetrances for these six genotypes are always equal to zero regardless of the heritability. On the other hand, genotypes AAbb, AaBb and aaBB have nonzero penetrances (see Methods for details). Therefore, different heritability levels certainly lead to different penetrances for genotypes AAbb, AaBb and aaBB. However, the ratios between the penetrances of these three genotypes are fixed and independent of the heritability. This model description can be generalised to cover the other multilocus interaction models. In addition, the maximum penetrance in any twolocus or multilocus interaction models always stays below 0.1 even though the heritability is at the highest level (see Methods for details). This means that case samples are always oversampled from affected individuals to achieve a balanced casecontrol data set. Since all explored heritability levels lead to the same case oversampling pattern, the simulated data sets of which the only primary difference being the heritability levels are indistinguishable from one another. This leads to the result similarities in interaction problems with the same number of SNPs in the data set, sample size and number of causative SNPs but different levels of heritability as shown in Figures 3, 4, 5, 6, 7 and 8. The result trend is also independent of the number of simulated data sets used in the benchmark trial.
In a permutation test, the ability to differentiate between two pvalues is influenced by the number of permutation replicates. With t permutation replicates, the test declares an actual pvalue that is less than 1/t to be zero. During the progressive search for the best ensemble, the inclusion of a new twoSNP unit is accepted if this inclusion does not worsen the current result. If the number of permutation replicates is too low, the search may include erroneous twoSNP units that are irrelevant to the correct association explanation. The analysis is confirmed as the number of output SNPs from 2LOmb is equal to the number of causative SNPs in most of simulation results. This phenomenon suggests that the number of permutation replicates employed in this study (t = 10,000) is high enough to screen off most of the erroneous twoSNP units. In other words, the inclusion of these erroneous twoSNP units leads to an increase in the pvalue by at least 1/t. Nonetheless, the fact that 2LOmb results are not entirely free from erroneous SNPs suggests that there are erroneous twoSNP units with extremely small pvalues. It is advisable to perform a genotype relative risk calculation for the elimination of erroneous SNPs. If the presence of an erroneous twoSNP unit is suspected, its result on twolocus genotype relative risk would not be as significant as that from the other twoSNP units in the ensemble. Alternatively, an additional means for further SNP screening by other techniques such as MDR is also recommended. The chance of erroneous SNP discovery would be further minimised by employing two consecutive attribute selection techniques. The same concept has been adopted for the implementation of MDR software, in which many additional filters including a χ^{2} test, an odds ratio test, ReliefF and TuRF are available for SNP screening prior to the MDR analysis.
Six genotypes in Figure 10 namely AABB, AABb, AaBB, Aabb, aaBb and aabb are protective genotypes. In other words, a sample with one of these six genotypes is a control sample because the penetrances for these genotypes are zero. It is also noted that the control samples with all nine genotypes precisely follow the distribution as jointly described by independent singlelocus genotype distribution from loci A and B. In contrast, three remaining genotypes in Figure 10 namely AAbb, AaBb and aaBB are labelled as diseasepredisposing genotypes because the majority of samples with these three genotypes are case samples. Samples with these genotypes may be either case or control samples because the penetrances for these genotypes are between zero and one. In fact, the probabilities that persons with these genotypes to have the disease are quite low since the penetrances for these genotypes are small. However, case samples must be oversampled from affected individuals to ensure a balanced casecontrol data set because the disease prevalence for this twolocus interaction model is only 0.004975. In addition, each case sample must contain one of these three genotypes because the penetrances for the other genotypes are zero. As a result, the case samples with these genotypes do not follow the same twolocus genotype distribution as in the control samples. With six genotypes being exclusively specific to control samples and the majority of three remaining genotypes being found in case samples, the MDR prediction accuracy for the twolocus interaction model is high. This explanation can also be generalised to cover the MDR results from the other multilocus interaction data sets.
Computational time required by 2LOmb, a combined 2LOmb and MDR approach, and direct MDR analysis to detect interactions in simulated data sets with different sizes and different numbers of causative SNPs.
Computational Time Required by Each Approach (sec)  

2LOmb  2LOmb+MDR  MDR  
Number of Causative SNPs  Sample Size  1,000 SNPs  2,000 SNPs  4,000 SNPs  1,000 SNPs  2,000 SNPs  4,000 SNPs  100 SNPs 
2  400  15  37  135  17  39  137  7,656 
800  21  59  224  23  61  226  15,990  
1,600  36  106  400  38  108  402  31,222  
3  400  22  43  140  24  45  142  7,721 
800  30  65  229  32  67  231  16,206  
1,600  50  115  406  52  117  408  31,232  
4  400  32  55  150  34  57  152  7,841 
800  46  80  236  48  82  238  16,285  
1,600  70  133  419  72  135  421  31,637 
The simulated multilocus interaction problems in this article are based on the pure epistasis model by Culverhouse et al. [42]. It is possible to capture a number of multilocus interactions with marginal twolocus effects via a combination of twolocus analyses. However, there are many multilocus interaction scenarios without marginal twolocus effects. In such cases, 2LOmb and the exhaustive twolocus analysis technique are unable to detect interactions. Among the explored techniques, TuRF and MDR have a better chance of detection. Nonetheless, TuRF functions well only when the total number of SNPs in data is small and the sample size is large enough while the total number of SNPs in data affects the practicality of direct MDR analysis.
Every attribute selection technique has a limitation in terms of the maximum numbers of samples and attributes that it can handle. Singlelocus analysis techniques always have a higher limit than multilocus analysis techniques. Because attribute subset evaluation is usually integrated into multilocus analysis techniques, consequently the number of possible attribute subsets that can be explored is extremely large when the candidate attribute set is large. Together with a potentially large sample size, a higher computational requirement for multilocus analysis techniques is inevitable. As a result, the direct application of multilocus analysis techniques to a much larger data set than those presented in this article, which is usually considered in genomewide association studies, would be impractical. However, it is reasonable to expect that both marginal singlelocus and epistatic effects are present in any genomewide data sets. A multistage strategy that incorporates multiple techniques, designed for different detection modes, would be more suitable to handle large data. For instance, the marginal singlelocus effects should be the first priority and, as such, be detected by singlelocus analysis. Then, a special case of pure epistasis [2] or semipurely epistatic events, in which a SNP displaying a marginal singlelocus effect interacts with a SNP that exhibits no marginal singlelocus effect, should be considered. Many twolocus analysis techniques have been proven to be well suited to this type of epistasis [40, 50, 51]. Finally, the detection of pure epistasis is carried out in the last stage. With the reduction of SNPs from the first stage, the chance that some multilocus analysis techniques are applicable to the remaining SNPs increases. In addition to the multistage approach, a prior knowledge regarding the previously reported association can be exploited to select candidate genes based upon ontology and pathways. This practice is due to the necessity for the derivation of plausible interpretation. The screening for SNPs within or near candidate genes before the association detection also increases the chance that multilocus analysis techniques can be applied to the remaining data.
Testing with real data
2LOmb has been applied to study a type 2 diabetes mellitus (T2D) data set, collected and investigated by the Wellcome Trust Case Control Consortium (WTCCC) [3]. The data set consists of 1,999 case samples from affected individuals in the UK and 3,004 control samples, which are the results of a merging between 1,500 samples from the UK blood services and 1,504 samples from the 1958 British birth cohort. The original genomewide data set contains 500,568 SNPs that are obtained through the Affymetrix GeneChip 500 K Mapping Array Set. The SNP set is primarily reduced by screening for SNPs within and near 372 candidate genes collected by the Human Genome Epidemiology Network (HuGENet) [62]. These candidate genes cover genes from both positive and negative genetic association reports, in which studies are conducted in various ethnic groups and populations. The SNP set is further reduced by removing SNPs that exhibit strong evidence of genetic association via singlelocus analysis. The final SNP set contains 7,065 SNPs from 370 candidate genes. All SNPs in the reduced data set exhibit no marginal singlelocus effects (Bonferronicorrected χ^{2}'s pvalue > 0.05). Detailed description of the final SNP set is given in the supplement (see Additional file 1).
2LOmb identifies 11 intronic SNPs, which are located in four genes, from the reduced T2D data.
Gene  Chromosome and Location  TwoSNP Unit in the Ensemble 

PGM1 (phosphoglucomutase 1)  1p31  (rs2269241, rs3790857) 
(rs2269239, rs3790857)  
(rs3790857, rs2269238)  
LMX1A (LIM homeobox transcription factor 1, alpha)  1q22q23  (rs2348250, rs6702087) 
PARK2 (Parkinson disease (autosomal recessive, juvenile) 2, parkin)  6q25.2q27  (rs1893551, rs6924502) 
GYS2 (glycogen synthase 2 (liver))  12p12.2  (rs6487236, rs1871142) 
(rs1871142, rs10770836) 
Genotype relative risk evaluated from genotype distribution of SNPs in PGM1.
Frequency  

SNP  Genotype  Case  Control  Relative Risk  95% CI 
SNP1  0  0.6513  0.6528  0.9977  (0.95731.0399) 
1  0.3082  0.3076  1.0018  (0.92041.0905)  
2  0.0405  0.0396  1.0229  (0.77571.3488)  
SNP2  0  0.6493  0.6521  0.9957  (0.95521.0379) 
1  0.3087  0.3079  1.0024  (0.92091.0910)  
2  0.0420  0.0399  1.0519  (0.80061.3822)  
SNP3  0  0.6153  0.6568  0.9368  (0.89720.9782) 
1  0.3472  0.3056  1.1361  (1.04791.2316)  
2  0.0375  0.0376  0.9974  (0.74901.3281)  
SNP4  0  0.6638  0.6668  0.9956  (0.95641.0364) 
1  0.2991  0.2969  1.0074  (0.92371.0988)  
2  0.0370  0.0363  1.0202  (0.76361.3631)  
(SNP1, SNP3)  0 0  0.6043  0.6448  0.9372  (0.89660.9796) 
0 1  0.0470  0.0080  5.8858  (3.77309.1816)  
1 0  0.0110  0.0120  0.9183  (0.54201.5561)  
1 1  0.2961  0.2943  1.0064  (0.92221.0983)  
1 2  0.0010  0.0013  0.7514  (0.13784.0984)  
2 1  0.0040  0.0033  1.2022  (0.47533.0408)  
2 2  0.0365  0.0363  1.0064  (0.75231.3463)  
(SNP2, SNP3)  00  0.6038  0.6448  0.9364  (0.89580.9789) 
01  0.0455  0.0073  6.2159  (3.91549.8681)  
10  0.0115  0.0117  0.9875  (0.58531.6661)  
11  0.2966  0.2949  1.0058  (0.92181.0975)  
12  0.0005  0.0013  0.3757  (0.04203.3588)  
21  0.0050  0.0033  1.5028  (0.62663.6038)  
22  0.0370  0.0363  1.0202  (0.76361.3631)  
(SNP3, SNP4)  00  0.6138  0.6551  0.9369  (0.89710.9785) 
01  0.0015  0.0017  0.9017  (0.21573.7686)  
10  0.0500  0.0117  4.2936  (2.93406.2831)  
11  0.2971  0.2936  1.0121  (0.92741.1044)  
21  0.0005  0.0017  0.3006  (0.03512.5706)  
22  0.0370  0.0360  1.0297  (0.77021.3765) 
Haplotype relative risk evaluated from genotype distribution of SNPs in PGM1.
Frequency  

SNP  Allele and Haplotype  Case  Control  Relative Risk  95% CI 
SNP1  0  0.8054  0.8066  0.9985  (0.97911.0183) 
1  0.1946  0.1934  1.0061  (0.92741.0916)  
SNP2  0  0.8037  0.8061  0.9970  (0.97751.0168) 
1  0.1963  0.1939  1.0126  (0.93361.0982)  
SNP3  0  0.7889  0.8096  0.9744  (0.95500.9943) 
1  0.2111  0.1904  1.1087  (1.02401.2003)  
SNP4  0  0.8134  0.8152  0.9977  (0.97891.0170) 
1  0.1866  0.1848  1.0100  (0.92881.0982)  
(SNP1, SNP3)  0 0  0.7812  0.8019  0.9742  (0.95430.9945) 
0 1  0.0242  0.0047  5.1533  (3.39467.8231)  
1 0  0.0077  0.0077  1.0000  (0.63491.5752)  
1 1  0.1869  0.1857  1.0064  (0.92571.0941)  
(SNP2, SNP3)  00  0.7804  0.8017  0.9734  (0.95350.9938) 
01  0.0232  0.0044  5.3216  (3.45578.1949)  
10  0.0085  0.0079  1.0758  (0.69301.6701)  
11  0.1879  0.1861  1.0099  (0.92911.0977)  
(SNP3, SNP4)  00  0.7881  0.8086  0.9747  (0.95520.9946) 
01  0.0008  0.0010  0.7661  (0.19433.0205)  
10  0.0253  0.0067  3.7936  (2.63675.4582)  
11  0.1858  0.1837  1.0113  (0.92981.0999) 
Genotype relative risk evaluated from genotype distribution of SNPs in LMX1A.
Frequency  

SNP  Genotype  Case  Control  Relative Risk  95% CI 
SNP5  0  0.8429  0.8642  0.9754  (0.95260.9987) 
1  0.1531  0.1315  1.1642  (1.01401.3366)  
2  0.0040  0.0043  0.9248  (0.38402.2271)  
SNP6  0  0.8799  0.8492  1.0362  (1.01351.0594) 
1  0.1161  0.1465  0.7924  (0.68290.9193)  
2  0.0040  0.0043  0.9248  (0.38402.2271)  
(SNP5, SNP6)  00  0.8329  0.8299  1.0036  (0.97841.0295) 
01  0.0100  0.0343  0.2918  (0.18140.4695)  
10  0.0470  0.0193  2.4355  (1.76443.3618)  
11  0.1061  0.1119  0.9482  (0.80611.1153)  
22  0.0040  0.0040  1.0018  (0.41032.4465) 
Haplotype relative risk evaluated from genotype distribution of SNPs in LMX1A.
Frequency  

SNP  Allele and Haplotype  Case  Control  Relative Risk  95% CI 
SNP5  0  0.9195  0.9299  0.9887  (0.97741.0002) 
1  0.0805  0.0701  1.1494  (0.99971.3214)  
SNP6  0  0.9380  0.9224  1.0168  (1.00591.0279) 
1  0.0620  0.0776  0.7997  (0.68920.9280)  
(SNP5, SNP6)  00  0.9143  0.9124  1.0021  (0.98981.0145) 
01  0.0051  0.0175  0.2931  (0.18290.4697)  
10  0.0236  0.0100  2.3640  (1.71503.2585)  
11  0.0569  0.0601  0.9472  (0.80641.1127) 
Genotype relative risk evaluated from genotype distribution of SNPs in PARK2.
Frequency  

SNP  Genotype  Case  Control  Relative Risk  95% CI 
SNP7  0  0.4802  0.5110  0.9398  (0.88730.9954) 
1  0.4322  0.4041  1.0695  (1.00081.1429)  
2  0.0875  0.0849  1.0313  (0.85811.2395)  
SNP8  0  0.4892  0.4923  0.9937  (0.93801.0527) 
1  0.4087  0.4121  0.9917  (0.92671.0613)  
2  0.1021  0.0955  1.0682  (0.90091.2665)  
(SNP7, SNP8)  00  0.4492  0.4844  0.9275  (0.87260.9858) 
01  0.0285  0.0260  1.0982  (0.78411.5380)  
02  0.0025  0.0007  3.7569  (0.729619.3450)  
10  0.0400  0.0080  5.0092  (3.18567.8767)  
11  0.3792  0.3848  0.9854  (0.91691.0590)  
12  0.0130  0.0113  1.1492  (0.69181.9089)  
21  0.0010  0.0013  0.7514  (0.13784.0984)  
22  0.0865  0.0836  1.0358  (0.86061.2465) 
Haplotype relative risk evaluated from genotype distribution of SNPs in PARK2.
Frequency  

SNP  Allele and Haplotype  Case  Control  Relative Risk  95% CI 
SNP7  0  0.6963  0.7130  0.9766  (0.95151.0023) 
1  0.3037  0.2870  1.0582  (0.99501.1254)  
SNP8  0  0.6936  0.6984  0.9931  (0.96721.0198) 
1  0.3064  0.3016  1.0159  (0.95631.0793)  
(SNP7, SNP8)  00  0.6726  0.6937  0.9696  (0.94340.9966) 
01  0.0238  0.0194  1.2248  (0.93691.6011)  
10  0.0210  0.0048  4.4215  (2.89726.7480)  
11  0.2826  0.2822  1.0016  (0.93971.0675) 
Genotype relative risk evaluated from genotype distribution of SNPs in GYS2.
Frequency  

SNP  Genotype  Case  Control  Relative Risk  95% CI 
SNP9  0  0.6473  0.6575  0.9846  (0.94471.0262) 
1  0.3167  0.3046  1.0396  (0.95581.1308)  
2  0.0360  0.0379  0.9491  (0.71051.2679)  
SNP10  0  0.5933  0.6441  0.9211  (0.88050.9634) 
1  0.3712  0.3169  1.1713  (1.08391.2657)  
2  0.0355  0.0389  0.9119  (0.68281.2180)  
SNP11  0  0.6058  0.6142  0.9864  (0.94271.0321) 
1  0.3507  0.3352  1.0461  (0.96751.1310)  
2  0.0435  0.0506  0.8601  (0.66501.1126)  
(SNP9, SNP10)  00  0.5863  0.6335  0.9255  (0.88410.9689) 
01  0.0610  0.0240  2.5463  (1.91353.3885)  
10  0.0055  0.0107  0.5166  (0.26101.0224)  
11  0.3092  0.2919  1.0590  (0.97171.1540)  
12  0.0020  0.0020  1.0018  (0.28313.5456)  
21  0.0010  0.0010  1.0018  (0.16765.9902)  
22  0.0335  0.0370  0.9071  (0.67341.2218)  
(SNP10, SNP11)  00  0.5463  0.5922  0.9224  (0.87760.9695) 
01  0.0455  0.0506  0.8997  (0.69821.1593)  
02  0.0015  0.0013  1.1271  (0.25255.0304)  
10  0.0595  0.0220  2.7095  (2.01643.6408)  
11  0.3032  0.2823  1.0739  (0.98391.1722)  
12  0.0085  0.0126  0.6723  (0.38051.1877)  
21  0.0020  0.0023  0.8587  (0.25172.9296)  
22  0.0335  0.0366  0.9153  (0.67911.2336) 
Haplotype relative risk evaluated from genotype distribution of SNPs in GYS2.
Frequency  

SNP  Allele and Haplotype  Case  Control  Relative Risk  95% CI 
SNP9  0  0.8057  0.8098  0.9949  (0.97571.0146) 
1  0.1943  0.1902  1.0216  (0.94121.1087)  
SNP10  0  0.7789  0.8026  0.9705  (0.95050.9908) 
1  0.2211  0.1974  1.1201  (1.03671.2102)  
SNP11  0  0.7811  0.7818  0.9992  (0.97821.0205) 
1  0.2189  0.2182  1.0030  (0.92991.0818)  
(SNP9, SNP10)  00  0.7740  0.7967  0.9715  (0.95120.9922) 
01  0.0317  0.0131  2.4258  (1.83573.2056)  
10  0.0049  0.0059  0.8330  (0.48071.4433)  
11  0.1894  0.1843  1.0276  (0.94551.1169)  
(SNP10, SNP11)  00  0.7494  0.7692  0.9742  (0.95240.9965) 
01  0.0295  0.0334  0.8843  (0.70691.1061)  
10  0.0318  0.0126  2.5275  (1.90643.3511)  
11  0.1894  0.1848  1.0244  (0.94261.1134) 
Prediction accuracy of the best MDR model constructed from the 2LOmb output.
Description  Value 

SNP and Gene  rs2269241 (PGM1), rs3790857 (PGM1), rs1893551 (PARK2), rs6924502 (PARK2), rs1871142 (GYS2), rs10770836 (GYS2) 
Classification (Training) Accuracy  0.5709 
Prediction Accuracy  0.5402 
CrossValidation Consistency (CVC)  9/10 
Summary of prediction accuracy by MDR from early genetic association studies of T2D in a Korean population, a Han Chinese population from Taiwan, a female population from the US, and that from an early genetic association study of metabolic syndrome in an Italian population from the Centre East Coast Italy.
The four genes selected by 2LOmb regulate many pathways that involve in the disease development [70–72]. The genetic association studies involving these genes have been previously conducted in different populations. For instance, LMX1A has been chosen as a positional and biological candidate gene for a casecontrol study of T2D in Pima Indians [73]. This gene is chosen as a candidate because a linkage of T2D to chromosome 1q21q23 has been previously reported [74]. In addition, LMX1A is one of LIM homeobox genes that are expressed in pancreas and has been shown to activate insulin gene transcription. Although SNPs have been carefully selected from the entire gene, no association between these SNPs in LMX1A and T2D has been found in this ethnic group.
PARK2 is another candidate gene that is also selected for casecontrol studies, based on evidence from genomewide linkage analysis [75]. A linkage of T2D in an African American population to chromosome 6q24q27 has been previously identified [76]. Although PARK2 mainly involves in the development of Parkinson's disease, singlelocus analysis reveals strong evidence of association between SNPs, which are in the vicinity of SNPs identified by 2LOmb, and T2D in African Americans.
In contrast to LMX1A and PARK2, which are candidate genes in typical T2D casecontrol studies, GYS2 is considered in a study to identify genes responsible for troglitazoneassociated hepatotoxicity in Japaneses with T2D [77]. In other words, both case and control samples in the study are drawn from troglitazonetreated T2D patients, in which case patients exhibit an abnormal increase in alanine transaminase (ALT) and aspartate transaminase (AST) levels. GYS2 regulates starch and sucrose metabolism and an insulin signalling pathway. The selected SNPs in GYS2 are not found to associate with troglitazoneinduced hepatotoxicity.
Similar to the study of GYS2, the association study involving PGM1 is not carried out as a typical T2D casecontrol study. In fact, an attempt to identify association between PGM1 polymorphisms and obesity has been conducted among T2D affected individuals in Italy [78]. PGM1 regulates glycolysis and gluconeogenesis, starch and sucrose metabolism, galactose metabolism, a pentose phosphate pathway, and streptomycin biosynthesis. Isozyme polymorphisms [79, 80], which are defined through structural differences in PGM1 protein, are used instead of SNPs in the study where positive association is identified.
In summary, positive association has been reported from previous studies involving PARK2 in African Americans and PGM1 in Italians. In contrast, negative association has been reported from previous studies about LMX1A in Pima Indians and GYS2 in Japaneses. Both GYS2 and PGM1 regulate starch and sucrose metabolism while LMX1A and PARK2 govern insulin gene transcription and Parkinson's disease development, respectively. The above discussion strengthens the importance of conducting largescale association studies due to two main reasons. Firstly, a gene that does not contribute to the aetiology of a complex disease in one population may be important for association explanation in another population. Secondly, the absence of interacting candidate genes from a study may lead to negative association due to a lack of necessary genetic information. A twolocus interaction can occur between SNPs from genes that regulate one specific pathway [44] or between SNPs from genes that regulate different pathways [45]. Furthermore, a multilocus interaction may involve both SNPs from genes that regulate the same pathway and SNPs from genes that govern different pathways. Hence, candidate genes should be selected by considering all pathways that directly and indirectly contribute to the disease development.
This study produces evidence of association between 11 intronic SNPs in PGM1, LMX1A, PARK2 and GYS2, and T2D in a UK population. Although there are other independent genomewide T2D data sets, the association detection within these data using a similar methodology to the presented method has never been attempted because the methodology employed in the majority of genomewide association studies is based on singlelocus analysis [3, 81]. It is recalled that each SNP explored in the reduced T2D data set exhibits no marginal singlelocus effect. Hence, the most logical approach to confirm the possibility of replicating association results from the current study is to perform the same detection method on these independent data sets. This is certainly important to gain further understanding of the genetic role in T2D susceptibility.
Implementation
2LOmb is implemented in a C programming language. All functions within the program are written by the first author except the χ^{2} distribution function, which is taken from the Numerical Recipes in C [82]. The program can be compiled by Microsoft Visual Studio and GNU C compilers. The program has been successfully tested for the execution under Windows and Linux operating systems. The time required by 2LOmb to complete a problem containing n attributes is T (n) = = n!/((n  2)!2!) = n(n  1)/ 2. 2LOmb thus has the order of n^{2} time complexity (T (n) ∈ O(n^{2})). Consequently, 2LOmb can tackle problems in quadratic time. 2LOmb in its present form occupies one processor during the program execution. A parallel version of 2LOmb for genomewide data is under development. All results included in the study are collected from the execution of computer programs in a Beowulf cluster. The computational platform consists of 12 nodes. Each node is equipped with dual Xeon 2.8 GHz processors and 4GB of main memory. The Rocks Cluster Distribution is installed on all nodes.
Conclusion
In this article, a method for detecting epistatic multilocus interactions in casecontrol data is presented. The study focuses on pure epistasis [2], which cannot be detected via singlelocus analysis [42]. To overcome this difficulty, the proposed method performs an omnibus permutation test [54] on ensembles of twolocus analyses and is thus referred to as 2LOmb. The detection performance of 2LOmb is evaluated using both simulated and real data. From the simulation, 2LOmb produces a low falsepositive error when the tests on null data of no association are performed. Furthermore, 2LOmb can identify all causative SNPs and outperforms a simple exhaustive twolocus analysis technique, a set association approach (SAA) [11], a correlationbased feature selection (CFS) technique [14] and a tuned ReliefF (TuRF) technique [16] in various interaction scenarios with marginal twolocus effects. These scenarios are set up by varying the number of causative SNPs, the number of SNPs in data, the sample size and the heritability. ANOVA reveals that the number of SNPs in data and the sample size influence the number of erroneous SNPs appended to the correctlyidentified causative SNPs in the 2LOmb output. In contrast, the results from 2LOmb appear to be insensitive to the variation in heritability. After subjecting the data sets containing only SNPs that are screened by 2LOmb to multifactor dimensionality reduction (MDR) analysis [19], all erroneous SNPs are successfully removed. In addition, an insight into the MDR models is provided. 2LOmb is subsequently applied to a real casecontrol type 2 diabetes mellitus (T2D) data set, which is collected from a UK population by the Wellcome Trust Case Control Consortium (WTCCC) [3]. The original genomewide data set is first reduced by selecting only SNPs that locate within or near 372 candidate genes reported by the Human Genome Epidemiology Network (HuGENet) [62]. In addition, the selected SNPs must exhibit no marginal singlelocus effects. The final data set, which consists of 1,999 case samples and 3,004 control samples, contains 7,065 SNPs from 370 candidate genes. 2LOmb identifies 11 intronic SNPs that are associated with the disease. These SNPs are located in PGM1, LMX1A, PARK2 and GYS2. The 2LOmb result suggests that there is no interaction between each pair of the identified genes that can be described by purely epistatic twolocus interaction models. Moreover, there are no interactions between these four genes that can be described by purely epistatic multilocus interaction models with marginal twolocus effects. This evidence of genetic association for these four genes leads to an alternative explanation for the aetiology of T2D in the UK population. It also implies that SNPs from genomewide data which are usually discarded after singlelocus analysis confirms the null hypothesis of no association can still be useful for genetic association studies of complex diseases.
Methods
Pure epistasis model
Penetrances for a twolocus interaction model.
Penetrance of Genotype  

Genotype  BB  Bb  bb  Marginal Penetrance 
AA  f _{00}  f _{01}  f _{02}  M _{A 0} 
Aa  f _{10}  f _{11}  f _{12}  M _{A 1} 
aa  f _{20}  f _{21}  f _{22}  M _{A 2} 
Marginal Penetrance  M _{B 0}  M _{B 1}  M _{B 2}  K 
The search for feasible penetrance f_{ ij }that also maximises the heritability or other variancebased objectives can be treated as a constraint optimisation problem. Many algorithms including a double description method [42] and a genetic algorithm [83] have been proven to be suitable for the task.
Twolocus penetrances that lead to the maximum heritability (K) = 2K/(1  K) for K ∈ (0, 1/4].
Penetrance of Genotype  

Genotype  BB  Bb  bb 
AA  0  0  4K 
Aa  0  2K  0 
aa  4K  0  0 
Threelocus penetrances that lead to the maximum heritability (K) = 9K/(1  K) for K ∈ (0, 1/16].
Penetrance of Genotype  

Genotype  CC  Cc  cc  
BB  Bb  bb  BB  Bb  bb  BB  Bb  bb  
AA  0  0  16K  0  0  0  0  0  0 
Aa  0  0  0  0  4K  0  0  0  0 
aa  0  0  0  0  0  0  16K  0  0 
Fourlocus penetrances that lead to the maximum heritability (K) = 35K/(1  K) for K ∈ (0, 1/64].
Penetrance of Genotype  

Genotype  CC  Cc  cc  
DD  Dd  dd  DD  Dd  dd  DD  Dd  Dd  
BB  0  0  0  0  0  0  0  0  0  
AA  Bb  0  0  0  0  0  0  0  0  0 
bb  0  0  64K  0  0  0  0  0  0  
BB  0  0  0  0  0  0  0  0  0  
Aa  Bb  0  0  0  0  8K  0  0  0  0 
bb  0  0  0  0  0  0  0  0  0  
BB  0  0  0  0  0  0  64K  0  0  
aa  Bb  0  0  0  0  0  0  0  0  0 
bb  0  0  0  0  0  0  0  0  0 
Disease prevalence that gives the target maximum heritability of 0.01, 0.025 and 0.05 for two, three and fourlocus interaction models.
Prevalence (K)  

Model  (K) = 0.01  (K) = 0.025  (K) = 0.05 
Twolocus  0.004975  0.012346  0.024390 
Threelocus  0.001110  0.002770  0.005525 
Fourlocus  0.000286  0.000714  0.001427 
genomeSIM
genomeSIM is a simulation package for generating casecontrol samples in largescale and genomewide association studies [60]. The package is capable of producing many realistic scenarios, which can be observed in a population and genetic samples, including linkage disequilibrium, phenocopy and genotyping errors. The case/control status of each sample is determined from the penetrancebased genetic models or interaction models. As a result, the package can accommodate many epistasis models including the one proposed by Culverhouse et al. [42]. A data set can be produced via two modes: a populationbased simulation and a probabilitybased simulation. In the populationbased simulation, an initial population is generated according to the predefined allele frequency of each SNP. Then further generations are created by crossing the genotype strings within successive generations until the specified number of generations is reached. The resulting data set contains a populationdependent case and control samples that follow a forwardtime simulation strategy. In contrast, genotype strings are incrementally generated without any string crossing for only one generation in the probabilitybased simulation. The creation of new strings is terminated only when the desired numbers of case and control samples are obtained. In this study, the probabilitybased simulation is used to produce all case and control samples where the simulation parameter setting is given in the supplement (see Additional file 2). genomeSIM is available upon request to Scott M. Dudek at the Vanderbilt University dudek@chgr.mc.vanderbilt.edu.
Set association approach
A set association approach (SAA) is an association detection technique based on an omnibus permutation test on sets of candidate SNPs [11]. The test captures information about genotyping errors, deviation from HardyWeinberg equilibrium (HWE) and allelic association. In the first step, the genotype distribution for each SNP in the control samples is checked for HWE. Then, the number of SNPs that is to be excluded from the study (n_{ d }) is set to the number of SNPs in the control samples that deviate from HWE. Two test statistics are subsequently calculated for each SNP: an allelic association statistic and a statistic for the deviation from HWE of each SNP in the case samples. The allelic association statistic is a χ^{2} statistic which is calculated from the contingency table of alleles or genotypes with disease status. On the other hand, a χ^{2} statistic for the deviation from HWE of each SNP in the case samples indicates the level of association. A large deviation from the equilibrium usually signifies strong association between a SNP and the disease. However, an excessively large deviation may be the result of genotyping errors. n_{ d }SNPs with largest test statistics for the deviation from HWE are hence excluded from the consideration.
The test statistics for the allelic association and deviation from HWE are multiplied together to form a single S statistic for each remaining SNP. SNPs are then ranked according to their S statistics. A preset number of SNPs with highest ranks are considered for association. The first candidate SNP set contains only the SNP with the highest rank (the highest S statistic). The pvalue for this first set is determined from a permutation simulation where the case and control labels are randomly permuted while the numbers of case and control samples remain unchanged. In each permutation replicate, a new genotype contingency table is constructed and a new S statistic is subsequently obtained. The pvalue is given by the fraction of permutation replicates with an S statistic greater than or equal to the S statistic from the original data. The second candidate SNP set consists of the first two SNPs in the rank list. The test statistic for this SNP set is the sum of S statistics from both SNPs. The pvalue for the second candidate SNP set is also obtained through the permutation simulation. By progressively adding the remaining SNP with the highest rank to the previously considered candidate set and performing the permutation simulation, pvalues for all candidate SNP sets are estimated. The sizes of candidate SNP sets have the range of one to the preset number. Among all candidate sets, the SNP set that best describes genetic association has the lowest pvalue.
Since multiple hypotheses are postulated during the construction of candidate SNP sets, the global pvalue for the selected candidate set must be evaluated. This is achieved through a permutation simulation in which the current raw pvalue for the chosen candidate set is now used as the test statistic. The existing permutation replicates, created for the early estimation of the raw pvalue, can be reused and a nested permutation simulation is hence avoided. In this study, the maximum allowable size of the candidate SNP set is the total number of available SNPs while the number of permutation replicates for pvalue evaluation is set to 10,000. The allelic association statistic employed in the study is the χ^{2} statistic that is obtained through the contingency table of genotypes with disease status. A PASCAL program for the set association approach can be obtained from the website for S statistic in gene mapping [84].
Correlationbased feature selection technique
where Merit_{ F }is the heuristic merit of an n_{ c }attribute subset F, _{ cf }is the average featureclass correlation and _{ ff }is the average featurefeature intercorrelation. An attribute subset receives a high merit score if it contains features that are highly correlated with the class and at the same time have low intercorrelation among one another. An application of a bestfirst search for the best subset identification is carried out to avoid searching through all possible attribute subsets. CFS has been integrated into a Weka package [85, 86].
Tuned ReliefF
A tuned ReliefF (TuRF) algorithm is a ranking algorithm for identifying genetic markers which are important in casecontrol classification [16]. TuRF is built on a ReliefF engine [15]. ReliefF randomly picks a sample from the (casecontrol) data and identifies its n_{ k }nearest neighbours from the same class and another n_{ k }nearest neighbours from the opposite class. The attribute valuesthe genotypes in this applicationof the neighbour samples are compared to that of the randomly picked sample and are subsequently used to update the relevance score for each attribute (genetic marker). This process is repeated for a specified number of samples, which is limited by the total sample size. The rationale of ReliefF is that an attribute which is important for the classification should have different values for samples from different classes and have the same value for samples from the same class. The relevance score of an attribute have a range from 1 (not relevant) to +1 (highly relevant). TuRF exploits the capability of ReliefF by repeatedly executing ReliefF and removing a portion of worst attributes at the end of each execution. This leads to the reevaluation of remaining attributes and, hence, reduces the effects of attribute noise on the attribute screening. In this study, the number of repetitions for random sample picking in the ReliefF part is equal to the total number of casecontrol samples while the neighbourhood size (n_{ k }) for the relevance score calculation is set to ten. Furthermore, the worst 1% of SNPs is removed at the end of each ReliefF iteration (TuRF 1%). TuRF has been integrated into the current distribution of multifactor dimensionality reduction (MDR) software.
Multifactor dimensionality reduction
The prediction accuracy of the decision table is subsequently evaluated by counting the numbers of case and control samples in the testing fold that their disease status can correctly be identified using the constructed decision rules. The process of decision table construction and evaluation must be cycled through all or some of possible  1 combinations where n_{ m }is the total number of available markers in the study. The best genetic marker combination is determined from two criteria: prediction accuracy and crossvalidation consistency. Each time that a testing fold is used for the prediction accuracy determination, the accuracy of the interesting marker combination model is compared with that from other models that also contain the same number of markers. The model that consistently ranks the first in comparison to other choices with the same number of markers has high crossvalidation consistency. Prediction accuracy is the main criterion for decision making while crossvalidation consistency is only used as an auxiliary measure. Crossvalidation consistency generally confirms that the high rank model can consistently be identified regardless of how the samples are divided for crossvalidation. In a situation where two or more models with different number of markers are equally good in terms of prediction accuracy and crossvalidation consistency, the most parsimonious modelthe combination with the least number of markersis chosen as the best model.
After the best model has been selected, a permutation test is used to assess the probability of obtaining prediction accuracy that is at least as large as or larger than that observed in the original data from randomised data. This represents the probability that the null hypothesis of no association is true. Each permutation replicate is constructed by randomly assigning the case/control status to each sample with the numbers of case and control samples remaining fixed. MDR analysis is subsequently carried out to obtain the prediction accuracy of each permutation replicate. The empirical pvalue is denoted by the fraction of permutation replicates with the prediction accuracy greater than or equal to the prediction accuracy obtained from the original data. MDR software, which incorporates many additional features including interaction visualisation via dendrograms and genetic marker screening via a χ^{2} test, an odds ratio test, ReliefF and TuRF, is available from its homepage [87].
JLIN
JLIN or a Java LINkage disequilibrium plotter is a computer program for visualisation of linkage disequilibrium analysis [63]. The program is capable of displaying many statistical measures including D'[64] and r^{2}[65]. The program is publicly available from the Centre for Genetic Epidemiology and Biostatistics, University of Western Australia [88].
Interaction dendrogram
An interaction dendrogram is a graphical tool for the visualisation of relationships among attributes (SNPs) [68, 69]. The interaction dendrogram is constructed via hierarchical clustering analysis and is embedded into MDR software [87]. The dendrogram illustrates the entropybased interaction between attributes by displaying interacting or related attributes closely together as adjacent leaves in a tree. At the same time, independent attributes are placed far apart from one another. In addition, the conclusion regarding whether the interaction between attributes is synergistic or redundancy is present can be deduced.
Availability and requirements
The 2LOmb program for Windows platforms and examples of simulated data are available at http://code.google.com/p/nachol/w/list.
Authors' information
WW is a Ph.D. student at the Department of Electrical Engineering, Faculty of Engineering, King Mongkut's University of Technology North Bangkok. He also received his B.Eng. and M.Eng. degrees in electrical engineering from King Mongkut's University of Technology North Bangkok. His current research interests include machine learning, evolutionary computation and bioinformatics.
AA is a Ph.D. student at the Department of Immunology, Faculty of Medicine Siriraj Hospital, Mahidol University. He also received his B.Sc. degree in pharmacy from Mahidol University. His current research interests include human genetics, genetic epidemiology, population genetics and bioinformatics.
TP is a Ph.D. student at the Department of Electrical Engineering, Faculty of Engineering, King Mongkut's University of Technology North Bangkok. He also received his B.Eng. and M.Eng. degrees in production engineering from King Mongkut's University of Technology North Bangkok. His current research interests include evolutionary multiobjective optimisation and machine learning.
SS is a parttime research assistant at the Department of Electrical Engineering, Faculty of Engineering, King Mongkut's University of Technology North Bangkok. He received his B.Eng. and M.Eng. degrees in electrical engineering from Thammasat University and King Mongkut's University of Technology North Bangkok, respectively. His current research interests include machine learning and genetic epidemiology.
CL is the Head of Division of Molecular Genetics at the Department of Research and Development, Faculty of Medicine Siriraj Hospital, Mahidol University. He also received his M.D. degree from Mahidol University. His current research interests include human genetics and genetic diseases.
NC is an associate professor of electrical engineering at King Mongkut's University of Technology North Bangkok and an adjunct professor of genetic epidemiology at Mahidol University. He received his B.Eng. and Ph.D. degrees from the Department of Automatic Control and Systems Engineering, University of Sheffield. His current research interests include evolutionary computation, machine learning and genetic epidemiology.
List of abbreviations
 2LOmb:

omnibus permutation test on ensembles of twolocus analyses
 ALT:

alanine transaminase
 ANOVA:

analysis of variance
 AST:

aspartate transaminase
 CFS:

correlationbased feature selection
 CI:

confidence interval
 CVC:

crossvalidation consistency
 DIO2 :

deiodinase, iodothyronine, type II
 E2LA:

exhaustive twolocus analysis
 EGFR :

epidermal growth factor receptor (erythroblastic leukemia viral (verbb) oncogene homolog, avian)
 FAMHAP:

software for singlemarker analysis and joint analysis of unphased genotype data from tightly linked markers (haplotype analysis)
 FUSION:

FinlandUnited States Investigation of NIDDM Genetics
 genomeSIM:

simulation package for generating casecontrol samples in largescale and genomewide association studies
 GYS2 :

glycogen synthase 2 (liver)
 HNF4A :

hepatocyte nuclear factor 4, alpha
 HuGENet:

Human Genome Epidemiology Network
 HWE:

HardyWeinberg equilibrium
 JLIN:

Java LINkage disequilibrium plotter
 KCNJ11 :

potassium inwardlyrectifying channel, subfamily J, member 11
 LD:

linkage disequilibrium
 LIM domains:

protein structural domains that are named after their initial discovery in the proteins Lin11, Isl1 and Mec3
 LMX1A :

LIM homeobox transcription factor 1, alpha
 MDR:

multifactor dimensionality reduction
 NIDDM:

noninsulindependent diabetes mellitus
 PARK2 :

Parkinson disease (autosomal recessive, juvenile) 2, parkin
 PGM1 :

phosphoglucomutase 1
 PPARG :

peroxisome proliferatoractivated receptor gamma
 RXRG :

retinoid X receptor, gamma
 SAA:

set association approach
 SNP:

single nucleotide polymorphism
 T2D:

type 2 diabetes mellitus
 TuRF:

tuned ReliefF
 UCP2 :

uncoupling protein 2 (mitochondrial, proton carrier)
 Weka:

Waikato environment for knowledge analysis
 WTCCC:

Wellcome Trust Case Control Consortium.
Declarations
Acknowledgements
The authors are extremely grateful to four anonymous reviewers and Dr. Danielle Talbot for their valuable comments and suggestions, which have contributed a lot towards improving the content and presentation of this article. The authors are also extremely grateful to Dr. Kuntinee Maneeratana, Dr. Pensiri Tongpadungrod, Dr. Graeme James Sheppard and Mr. Edward L. Karas for their efforts on proofreading the manuscript. The authors wish to thank Dr. Marong Phadoongsidhi and Dr. Vara Varavithya and for their assistance on processing the T2D data and providing an access to the Beowulf cluster, respectively. The authors also wish to thank Dr. Sissades Tongsima for his assistance on testing of multiple hypotheses, providing many suggestions for the manuscript revision and allowing the usage of computer facilities at the Genome Institute, National Center for Genetic Engineering and Biotechnology (BIOTEC), National Science and Technology Development Agency (NSTDA) during the algorithm development. The authors acknowledge Scott M. Dudek at the Vanderbilt University for providing an access to the genomeSIM package. Furthermore, the authors acknowledge Dr. Audrey Duncanson and the Wellcome Trust Case Control Consortium for granting an access to the genomewide casecontrol data. WW, AA and TP were supported by the Thailand Research Fund (TRF) through the Royal Golden Jubilee Ph.D. Programme (Grant No. PHD/1.E.KN.49/A.1, PHD/4.I.MU.45/C.1 and PHD/1.E.KN.50/A.1, respectively). CL was supported by the Mahidol Research Grant. NC was supported by the Thailand Research Fund and the National Science and Technology Development Agency through the NSTDA Chair Professor Grant.
Authors’ Affiliations
References
 Risch N, Merikangas K: The future of genetic studies of complex human diseases. Science 1996, 273: 1516–1517. 10.1126/science.273.5281.1516View ArticlePubMedGoogle Scholar
 Musani SK, Shriner D, Liu N, Feng R, Coffey CS, Yi N, Tiwari HK, Allison DB: Detection of gene × gene interactions in genomewide association studies of human population data. Hum Hered 2007, 63: 67–84. 10.1159/000099179View ArticlePubMedGoogle Scholar
 The Wellcome Trust Case Control Consortium: Genomewide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature 2007, 447: 661–678. 10.1038/nature05911PubMed CentralView ArticleGoogle Scholar
 The GAIN Collaborative Research Group: New models of collaboration in genomewide association studies: the Genetic Association Information Network. Nat Genet 2007, 39: 1045–1051. 10.1038/ng2127View ArticleGoogle Scholar
 Heidema AG, Boer JMA, Nagelkerke N, Mariman ECM, van der A DL, Feskens EJM: The challenge for genetic epidemiologists: how to analyze large numbers of SNPs in relation to complex diseases. BMC Genet 2006, 7: 23. 10.1186/14712156723PubMed CentralView ArticlePubMedGoogle Scholar
 Motsinger AA, Ritchie MD, Reif DM: Novel methods for detecting epistasis in pharmacogenomics studies. Pharmacogenomics 2007, 8: 1229–1241. 10.2217/14622416.8.9.1229View ArticlePubMedGoogle Scholar
 Moore JH, Gilbert JC, Tsai CT, Chiang FT, Holden T, Barney N, White BC: A flexible computational framework for detecting, characterizing, and interpreting statistical patterns of epistasis in genetic studies of human disease susceptibility. J Theor Biol 2006, 241: 252–261. 10.1016/j.jtbi.2005.11.036View ArticlePubMedGoogle Scholar
 Saeys Y, Inza I, Larrañaga P: A review of feature selection techniques in bioinformatics. Bioinformatics 2007, 23: 2507–2517. 10.1093/bioinformatics/btm344View ArticlePubMedGoogle Scholar
 Lewis CM: Genetic association studies: design, analysis and interpretation. Brief Bioinform 2002, 3: 146–153. 10.1093/bib/3.2.146View ArticlePubMedGoogle Scholar
 Montana G: Statistical methods in genetics. Brief Bioinform 2006, 7: 297–308. 10.1093/bib/bbl028View ArticlePubMedGoogle Scholar
 Hoh J, Wille A, Ott J: Trimming, weighting, and grouping SNPs in human casecontrol association studies. Genome Res 2001, 11: 2115–2119. 10.1101/gr.204001PubMed CentralView ArticlePubMedGoogle Scholar
 Potter DM: Omnibus permutation tests of the association of an ensemble of genetic markers with disease in casecontrol studies. Genet Epidemiol 2006, 30: 438–446. 10.1002/gepi.20155View ArticlePubMedGoogle Scholar
 Chapman J, Clayton D: Detecting association using epistatic information. Genet Epidemiol 2007, 31: 894–909. 10.1002/gepi.20250View ArticlePubMedGoogle Scholar
 Hall MA, Holmes G: Benchmarking attribute selection techniques for discrete class data mining. IEEE Trans Knowl Data Eng 2003, 15: 1437–1447. 10.1109/TKDE.2003.1245283View ArticleGoogle Scholar
 RobnikŠikonja M, Kononenko I: Theoretical and empirical analysis of ReliefF and RReliefF. Mach Learn 2003, 53: 23–69. 10.1023/A:1025667309714View ArticleGoogle Scholar
 Moore JH, White BC: Tuning ReliefF for genomewide genetic analysis. In Evolutionary Computation, Machine Learning and Data Mining in Bioinformatics. Edited by: Marchiori E, Moore JH, Rajapakse JC. Berlin, Heidelberg: Springer; 2007:166–175. [Goos G, Hartmanis J, van Leeuwen J (Founding and Former Series Editors): Lecture Notes in Computer Science, vol 4447]. [Goos G, Hartmanis J, van Leeuwen J (Founding and Former Series Editors): Lecture Notes in Computer Science, vol 4447].View ArticleGoogle Scholar
 Nelson MR, Kardia SLR, Ferrell RE, Sing CF: A combinatorial partitioning method to identify multilocus genotypic partitions that predict quantitative trait variation. Genome Res 2001, 11: 458–470. 10.1101/gr.172901PubMed CentralView ArticlePubMedGoogle Scholar
 Culverhouse R, Klein T, Shannon W: Detecting epistatic interactions contributing to quantitative traits. Genet Epidemiol 2004, 27: 141–152. 10.1002/gepi.20006View ArticlePubMedGoogle Scholar
 Ritchie MD, Hahn LW, Roodi N, Bailey LR, Dupont WD, Parl FF, Moore JH: Multifactordimensionality reduction reveals highorder interactions among estrogenmetabolism genes in sporadic breast cancer. Am J Hum Genet 2001, 69: 138–147. 10.1086/321276PubMed CentralView ArticlePubMedGoogle Scholar
 Hahn LW, Ritchie MD, Moore JH: Multifactor dimensionality reduction software for detecting genegene and geneenvironment interactions. Bioinformatics 2003, 19: 376–382. 10.1093/bioinformatics/btf869View ArticlePubMedGoogle Scholar
 Bush WS, Dudek SM, Ritchie MD: Parallel multifactor dimensionality reduction: a tool for the largescale analysis of genegene interactions. Bioinformatics 2006, 22: 2173–2174. 10.1093/bioinformatics/btl347View ArticlePubMedGoogle Scholar
 Chung Y, Lee SY, Elston RC, Park T: Odds ratio based multifactordimensionality reduction method for detecting genegene interactions. Bioinformatics 2007, 23: 71–76. 10.1093/bioinformatics/btl557View ArticlePubMedGoogle Scholar
 Bush WS, Edwards TL, Dudek SM, McKinney BA, Ritchie MD: Alternative contingency table measures improve the power and detection of multifactor dimensionality reduction. BMC Bioinformatics 2008, 9: 238. 10.1186/147121059238PubMed CentralView ArticlePubMedGoogle Scholar
 Lou XY, Chen GB, Yan L, Ma JZ, Mangold JE, Zhu J, Elston RC, Li MD: A combinatorial approach to detecting genegene and geneenvironment interactions in family studies. Am J Hum Genet 2008, 83: 457–467. 10.1016/j.ajhg.2008.09.001PubMed CentralView ArticlePubMedGoogle Scholar
 Edwards TL, Lewis K, Velez DR, Dudek SM, Ritchie MD: Exploring the performance of multifactor dimensionality reduction in large scale SNP studies and in the presence of genetic heterogeneity among epistatic disease models. Hum Hered 2009, 67: 183–192. 10.1159/000181157PubMed CentralView ArticlePubMedGoogle Scholar
 Mechanic LE, Luke BT, Goodman JE, Chanock SJ, Harris CC: Polymorphism Interaction Analysis (PIA): a method for investigating complex genegene interactions. BMC Bioinformatics 2008, 9: 146. 10.1186/147121059146PubMed CentralView ArticlePubMedGoogle Scholar
 Liang KH, Hwang Y, Shao WC, Chen EY: An algorithm for model construction and its applications to pharmacogenomic studies. J Hum Genet 2006, 51: 751–759. 10.1007/s1003800600162View ArticlePubMedGoogle Scholar
 EstradaGil JK, FernándezLópez JC, HernándezLemus E, SilvaZolezzi I, HidalgoMiranda A, JiménezSánchez G, VallejoClemente EE: GPDTI: a Genetic Programming Decision Tree Induction method to find epistatic effects in common complex diseases. Bioinformatics 2007, 23: i167i174. 10.1093/bioinformatics/btm205View ArticlePubMedGoogle Scholar
 Nunkesser R, Bernholt T, Schwender H, Ickstadt K, Wegener I: Detecting highorder interactions of single nucleotide polymorphisms using genetic programming. Bioinformatics 2007, 23: 3280–3288. 10.1093/bioinformatics/btm522View ArticlePubMedGoogle Scholar
 Lunetta KL, Hayward LB, Segal J, van Eerdewegh P: Screening largescale association study data: exploiting interactions using random forests. BMC Genet 2004, 5: 32. 10.1186/14712156532PubMed CentralView ArticlePubMedGoogle Scholar
 Bureau A, Dupuis J, Falls K, Lunetta KL, Hayward B, Keith TP, van Eerdewegh P: Identifying SNPs predictive of phenotype using random forests. Genet Epidemiol 2005, 28: 171–182. 10.1002/gepi.20041View ArticlePubMedGoogle Scholar
 Chen X, Liu CT, Zhang M, Zhang H: A forestbased approach to identifying gene and genegene interactions. Proc Natl Acad Sci USA 2007, 104: 19199–19203. 10.1073/pnas.0709868104PubMed CentralView ArticlePubMedGoogle Scholar
 Ritchie MD, White BC, Parker JS, Hahn LW, Moore JH: Optimization of neural network architecture using genetic programming improves detection and modeling of genegene interactions in studies of human diseases. BMC Bioinformatics 2003, 4: 28. 10.1186/14712105428PubMed CentralView ArticlePubMedGoogle Scholar
 MotsingerReif AA, Dudek SM, Hahn LW, Ritchie MD: Comparison of approaches for machinelearning optimization of neural networks for detecting genegene interactions in genetic epidemiology. Genet Epidemiol 2008, 32: 325–340. 10.1002/gepi.20307View ArticlePubMedGoogle Scholar
 Cordell HJ: Epistasis: what it means, what it doesn't mean, and statistical methods to detect it in humans. Hum Mol Genet 2002, 11: 2463–2468. 10.1093/hmg/11.20.2463View ArticlePubMedGoogle Scholar
 Wilson SR: Epistasis. In Nature Encyclopedia of the Human Genome. Volume 2. Edited by: Cooper DN. London: Nature Publishing Group; 2004:317–320.Google Scholar
 Neuman RJ, Rice JP: Twolocus models of disease. Genet Epidemiol 1992, 9: 347–365. 10.1002/gepi.1370090506View ArticlePubMedGoogle Scholar
 Schork NJ, Boehnke M, Terwilliger JD, Ott J: Twotraitlocus linkage analysis: a powerful strategy for mapping complex genetic traits. Am J Hum Genet 1993, 53: 1127–1136.PubMed CentralPubMedGoogle Scholar
 Li W, Reich J: A complete enumeration and classification of twolocus disease models. Hum Hered 2000, 50: 334–349. 10.1159/000022939View ArticlePubMedGoogle Scholar
 Marchini J, Donnelly P, Cardon LR: Genomewide strategies for detecting multiple loci that influence complex diseases. Nat Genet 2005, 37: 413–417. 10.1038/ng1537View ArticlePubMedGoogle Scholar
 Hallgrímsdóttir IB, Yuster DS: A complete classification of epistatic twolocus models. BMC Genet 2008, 9: 17. 10.1186/14712156917PubMed CentralView ArticlePubMedGoogle Scholar
 Culverhouse R, Suarez BK, Lin J, Reich T: A perspective on epistasis: limits of models displaying no main effect. Am J Hum Genet 2002, 70: 461–471. 10.1086/338759PubMed CentralView ArticlePubMedGoogle Scholar
 Cho YM, Ritchie MD, Moore JH, Park JY, Lee KU, Shin HD, Lee HK, Park KS: Multifactordimensionality reduction shows a twolocus interaction associated with type 2 diabetes mellitus. Diabetologia 2004, 47: 549–554. 10.1007/s0012500313213View ArticlePubMedGoogle Scholar
 Hsieh CH, Liang KH, Hung YJ, Huang LC, Pei D, Liao YT, Kuo SW, Bey MSJ, Chen JL, Chen EY: Analysis of epistasis for diabetic nephropathy among type 2 diabetic patients. Hum Mol Genet 2006, 15: 2701–2708. 10.1093/hmg/ddl203View ArticlePubMedGoogle Scholar
 Qi L, van Dam RM, Asselbergs FW, Hu FB: Genegene interactions between HNF4A and KCNJ11 in predicting type 2 diabetes in women. Diabet Med 2007, 24: 1187–1191. 10.1111/j.14645491.2007.02255.xView ArticlePubMedGoogle Scholar
 Zhang Z, Zhang S, Wong MY, Wareham NJ, Sha Q: An ensemble learning approach jointly modeling main and interaction effects in genetic association studies. Genet Epidemiol 2008, 32: 285–300. 10.1002/gepi.20304PubMed CentralView ArticlePubMedGoogle Scholar
 Fiorito M, Torrente I, De Cosmo S, Guida V, Colosimo A, Prudente S, Flex E, Menghini R, Miccoli R, Penno G, Pellegrini F, Tassi V, Federici M, Trischitta V, Dallapiccola B: Interaction of DIO2 T92A and PPARγ2 P12A polymorphisms in the modulation of metabolic syndrome. Obesity 2007, 15: 2889–2895. 10.1038/oby.2007.343View ArticlePubMedGoogle Scholar
 Albrechtsen A, Castella S, Andersen G, Hansen T, Pedersen O, Nielsen R: A Bayesian multilocus association method: allowing for higherorder interaction in association studies. Genetics 2007, 176: 1197–1208. 10.1534/genetics.107.071696PubMed CentralView ArticlePubMedGoogle Scholar
 Zhang Y, Liu JS: Bayesian inference of epistatic interactions in casecontrol studies. Nat Genet 2007, 39: 1167–1173. 10.1038/ng2110View ArticlePubMedGoogle Scholar
 Evans DM, Marchini J, Morris AP, Cardon LR: Twostage twolocus models in genomewide association. PLoS Genet 2006, 2: e157. 10.1371/journal.pgen.0020157PubMed CentralView ArticlePubMedGoogle Scholar
 Ionita I, Man M: Optimal twostage strategy for detecting interacting genes in complex diseases. BMC Genet 2006, 7: 39.PubMedGoogle Scholar
 Gayán J, GonzálezPérez A, Bermudo F, Sáez ME, Royo JL, Quintas A, Galan JJ, Morón FJ, RamirezLorca R, Real LM, Ruiz A: A method for detecting epistasis in genomewide studies using casecontrol multilocus association analysis. BMC Genomics 2008, 9: 360. 10.1186/147121649360PubMed CentralView ArticlePubMedGoogle Scholar
 Heidema AG, Feskens EJM, Doevendans PAFM, Ruven HJT, van Houwelingen HC, Mariman ECM, Boer JMA: Analysis of multiple SNPs in genetic association studies: comparison of three multilocus methods to prioritize and select SNPs. Genet Epidemiol 2007, 31: 910–921. 10.1002/gepi.20251View ArticlePubMedGoogle Scholar
 Pesarin F: Multivariate Permutation Tests with Applications to Biostatistics. Chichester: Wiley; 2001.Google Scholar
 Fisher RA: Statistical Methods for Research Workers. 4th edition. London: Oliver and Boyd; 1932.Google Scholar
 Westfall PH, Young SS: ResamplingBased Multiple Testing: Examples and Methods for pvalue Adjustment. New York: John Wiley and Sons; 1993.Google Scholar
 Becker T, Schumacher J, Cichon S, Baur MP, Knapp M: Haplotype interaction analysis of unlinked regions. Genet Epidemiol 2005, 29: 313–322. 10.1002/gepi.20096View ArticlePubMedGoogle Scholar
 Herold C, Becker T: Genetic association analysis with FAMHAP: a major program update. Bioinformatics 2009, 25: 134–136. 10.1093/bioinformatics/btn581View ArticlePubMedGoogle Scholar
 Hardy GH: Mendelian proportions in a mixed population. Science 1908, 28: 49–50. 10.1126/science.28.706.49View ArticlePubMedGoogle Scholar
 Dudek SM, Motsinger AA, Velez DR, Williams SM, Ritchie MD: Data simulation software for wholegenome association and other studies in human genetics. In Proceedings of the Pacific Symposium on Biocomputing 2006: 3–7 January 2006; Maui. Edited by: Altman RB, Dunker AK, Hunter L, Murray T, Klein TE. Singapore: World Scientific; 2006:499–510.View ArticleGoogle Scholar
 Guyon I, Elisseeff A: An introduction to variable and feature selection. J Mach Learn Res 2003, 3: 1157–1182. 10.1162/153244303322753616Google Scholar
 Yu W, Gwinn M, Clyne M, Yesupriya A, Khoury MJ: A navigator for human genome epidemiology. Nat Genet 2008, 40: 124–125. 10.1038/ng0208124View ArticlePubMedGoogle Scholar
 Carter KW, McCaskie PA, Palmer LJ: JLIN: a java based linkage disequilibrium plotter. BMC Bioinformatics 2006, 7: 60. 10.1186/14712105760PubMed CentralView ArticlePubMedGoogle Scholar
 Lewontin RC: The interaction of selection and linkage. I. general considerations; heterotic models. Genetics 1964, 49: 49–67.PubMed CentralPubMedGoogle Scholar
 Hill WG, Robertson A: Linkage disequilibrium in finite populations. Theor Appl Genet 1968, 38: 226–231. 10.1007/BF01245622View ArticlePubMedGoogle Scholar
 Excoffier L, Slatkin M: Maximumlikelihood estimation of molecular haplotype frequencies in a diploid population. Mol Biol Evol 1995, 12: 921–927.PubMedGoogle Scholar
 Epstein MP, Satten GA: Inference on haplotype effects in casecontrol studies using unphased genotype data. Am J Hum Genet 2003, 73: 1316–1329. 10.1086/380204PubMed CentralView ArticlePubMedGoogle Scholar
 Jakulin A, Bratko I, Smrke D, Demšar J, Zupan B: Attribute interactions in medical data analysis. In Artificial Intelligence in Medicine. Edited by: Dojat M, Keravnou E, Barahona P. Berlin, Heidelberg: Springer; 2003:229–238. [Carbonell JG, Siekmann J (Series Editors): Lecture Notes in Artificial Intelligence, vol 2780]. [Carbonell JG, Siekmann J (Series Editors): Lecture Notes in Artificial Intelligence, vol 2780].View ArticleGoogle Scholar
 Jakulin A, Bratko I: Analyzing attribute dependencies. In Knowledge Discovery in Databases: PKDD 2003. Edited by: Lavrač N, Gamberger D, Todorovski L, Blockeel H. Berlin, Heidelberg: Springer; 2003:229–240. [Carbonell JG, Siekmann J (Series Editors): Lecture Notes in Artificial Intelligence, vol 2838]. [Carbonell JG, Siekmann J (Series Editors): Lecture Notes in Artificial Intelligence, vol 2838].View ArticleGoogle Scholar
 Kanehisa M, Goto S: KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res 2000, 28: 27–30. 10.1093/nar/28.1.27PubMed CentralView ArticlePubMedGoogle Scholar
 Kanehisa M, Goto S, Hattori M, AokiKinoshita K, Itoh M, Kawashima S, Katayama T, Araki M, Hirakawa M: From genomics to chemical genomics: new developments in KEGG. Nucleic Acids Res 2006, 34: D354D357. 10.1093/nar/gkj102PubMed CentralView ArticlePubMedGoogle Scholar
 Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y: KEGG for linking genomes to life and the environment. Nucleic Acids Res 2008, 36: D480D484. 10.1093/nar/gkm882PubMed CentralView ArticlePubMedGoogle Scholar
 Thameem F, Wolford JK, Wang J, German MS, Bogardus C, Prochazka M: Cloning, expression and genomic structure of human LMX1A , and variant screening in Pima Indians. Gene 2002, 290: 217–225. 10.1016/S03781119(02)005826View ArticlePubMedGoogle Scholar
 Hanson RL, Ehm MG, Pettitt DJ, Prochazka M, Thompson DB, Timberlake D, Foroud T, Kobes S, Baier L, Burns DK, Almasy L, Blangero J, Garvey WT, Bennett PH, Knowler WC: An autosomal genomic scan for loci linked to type II diabetes mellitus and bodymass index in Pima Indians. Am J Hum Genet 1998, 63: 1130–1138. 10.1086/302061PubMed CentralView ArticlePubMedGoogle Scholar
 Leak TS, Mychaleckyj JC, Smith SG, Keene KL, Gordon CJ, Hicks PJ, Freedman BI, Bowden DW, Sale MM: Evaluation of a SNP map of 6q24–27 confirms diabetic nephropathy loci and identifies novel associations in type 2 diabetes patients with nephropathy from an AfricanAmerican population. Hum Genet 2008, 124: 63–71. 10.1007/s0043900805237PubMed CentralView ArticlePubMedGoogle Scholar
 Sale MM, Freedman BI, Langefeld CD, Williams AH, Hicks PJ, Colicigno CJ, Beck SR, Brown WM, Rich SS, Bowden DW: A genomewide scan for type 2 diabetes in AfricanAmerican families reveals evidence for a locus on chromosome 6q. Diabetes 2004, 53: 830–837. 10.2337/diabetes.53.3.830View ArticlePubMedGoogle Scholar
 Watanabe I, Tomita A, Shimizu M, Sugawara M, Yasumo H, Koishi R, Takahashi T, Miyoshi K, Nakamura K, Izumi T, Matsushita Y, Furukawa H, Haruyama H, Koga T: A study to survey susceptible genetic factors responsible for troglitazoneassociated hepatotoxicity in Japanese patients with type 2 diabetes mellitus. Clin Pharmacol Ther 2003, 73: 435–455. 10.1016/S00099236(03)000146View ArticlePubMedGoogle Scholar
 GloriaBottini F, Magrini A, Antonacci E, La Torre M, Di Renzo L, De Lorenzo A, Bergamaschi A, Bottini E: Phosphoglucomutase genetic polymorphism and body mass. Am J Med Sci 2007, 334: 421–425. 10.1097/MAJ.0b013e3180a5e934View ArticlePubMedGoogle Scholar
 Spencer N, Hopkinson DA, Harris H: Phosphoglucomutase polymorphism in man. Nature 1964, 204: 742–745. 10.1038/204742a0View ArticlePubMedGoogle Scholar
 March RE, Putt W, Hollyoake M, Ives JH, Lovegrove JU, Hopkinson DA, Edwards YH, Whitehouse DB: The classical human phosphoglucomutase (PGM1) isozyme polymorphism is generated by intragenic recombination. Proc Natl Acad Sci USA 1993, 90: 10730–10733. 10.1073/pnas.90.22.10730PubMed CentralView ArticlePubMedGoogle Scholar
 Zeggini E, Scott LJ, Saxena R, Voight BF: Metaanalysis of genomewide association data and largescale replication identifies additional susceptibility loci for type 2 diabetes. Nat Genet 2008, 40: 638–645. 10.1038/ng.120PubMed CentralView ArticlePubMedGoogle Scholar
 Press WH, Teukolsky SA, Vetterling WT, Flannery BP: Numerical Recipes in C: The Art of Scientific Computing. 2nd edition. Cambridge: Cambridge University Press; 1992.Google Scholar
 Moore JH, Hahn LW, Ritchie MD, Thornton TA, White BC: Routine discovery of complex genetic models using genetic algorithms. Appl Soft Comput 2004, 4: 79–86. 10.1016/j.asoc.2003.08.003PubMed CentralView ArticlePubMedGoogle Scholar
 S Statistic in Gene Mapping[http://www.genemapping.cn/sumstat.html]
 Witten IH, Frank E: Data Mining: Practical Machine Learning Tools and Techniques. 2nd edition. San Francisco: Morgan Kaufmann; 2005.Google Scholar
 Weka 3: Data Mining Software in Java[http://www.cs.waikato.ac.nz/ml/weka/]
 Multifactor Dimensionality Reduction[http://www.multifactordimensionalityreduction.org/]
 JLIN: A Java Based Linkage Disequilibrium Plotter[http://www.genepi.org.au/jlin.html]
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.