BACKGROUND
single snp association study is a popular method to detect genes that are susceptible to human diseases.
although candidate gene approach that uses prior knowledge about its function is an efficient procedure, it could overlook genes whose function is unknown or ambiguous.
gwas using whole genome snps is thus an attractive solution to this issue.
many procedures proposed for gwas so far are based on marginal association between each snp and phenotype.
however, it has been pointed out that most susceptible-snps identified often show low or moderate effect size, and hence may explain only a few percentage of the genetic variance  <cit> .
the fact that, for several complex diseases, the recurrence risk ratio decreases quicker than 1/ <dig> as the relatedness decreases implies the involvement of non-additive interactions in their etiology  <cit> .
in addition to substantial contributions of rarer variants, gene-gene and gene-environment interaction would be one of strong candidates that can explain the missing heritability as well  <cit> .
without investigating such interactions, we therefore may miss genuine disease-susceptible loci  <cit> .
an effective and accurate method to search gene-gene interactions will utilize immediately original snp-gwas data to decipher such missing genetic components of complex human diseases.
in this paper, we tackle a development of powerful method for the genome-wide gene-gene interaction analysis using snps.
marchinni et al.
<cit>  suggest that the use of arbitrary single locus-disease association model under which some interaction effect is present might prevent from finding susceptible interaction effect.
therefore exhaustive search for interactions are needed rather than primary screening by marginal effect.
although individual hypothesis tests for all snp-snp pairs as frequently used in snp-gwas might be the simplest approach, it is however faced with multiple testing issues caused by  <dig> - <dig> , <dig> snps.
specifically, there is a difficulty in setting genome-wide significance level using, e.g., bonferroni correction, fprp  <cit> , or bayes factor  <cit> , leading to prohibitively conservative result because they fail to successfully incorporate the correlation structure between each hypothesis.
for instance, the total number of hypothesises for second-order gene-gene interaction is about  <dig> – <dig> in standard gwas data.
then bonferroni-corrected significance level must be considerably small by multiplying the nominal significance level such as  <dig>  by the correction factor less than 10− <dig>  no efficient and universal multiple testing method to deal with such the huge set of hypotheses having complicated correlation structure is proposed so far.
although multiple logistic regression may be the method that can incorporate the correlation structure between hypothesises by putting them as predictor variables, the so-called large p small n situation , precludes simultaneous inclusion of all snp-snp interactions in the logistic regression model, namely, no unique solution is determined.
there are several software programs for gene-gene interaction analysis including multifactor dimensionality reduction   <cit> , fast-epistasis option in plink  <cit> , tuning relief  <cit> , random jungle  <cit> , beam  <cit> , and boost  <cit> .
cordell  <cit>  carried out an extensive comparative study for these methods.
she reported that mdr and beam have computational difficulty in analyzing current scale in gwas and that random jungle is applicable to small scale dataset but not for whole-genome dataset.
plink-fast-epistasis and boost are designed for exhaustive search.
since both methods are based on hypothesis testing principles, as stated above, exact determination of significance level for whole-genome data is difficult due to the numerous multiple hypotheses in addition to the complicated correlation structure between them.
hypothesis testing methods  <cit>  could be too conservative and lead to false negative result, because some complicated correlation structure would commonly appear when considering interaction pairs.
the boost has been shown to outperform the plink--fast-epistasis in most of the interaction models considered in  <cit> .
the boost primarily performs screening of possible hypotheses using kirkwood superposition approximation , then applies standard multiple testing on the basis of the initial number of snp-snp pairs to the likelihood ratio statistic under the saturated logistic regression model having the degrees of freedom  of  <dig>  however, the test statistic may fail to follow the chi-squared distribution of df  <dig> under the presence of sparse cells, which causes too small type  <dig> error rate  <cit> .
because the sparse cell issue is unavoidable in genotypic interaction data, some solutions in development of summary statistic for gene-gene interaction are strongly necessitated.
our method is free from the sparse cell issue owing to the proposed dummy coding strategy, which generates  <dig> by  <dig> contingency table.
in this paper we develop a new method based on contemporary ultrahigh-dimensional variable selection approach, designated the sure independence screening   <cit>  rather than hypothesis testing.
the principle of sis is based on the same as that of machine learning methods  <cit>  exploring the optimal model that can explain sufficiently the data in the most parsimonious way; its methodological feature to avoid over- or under-fitting enables us flexibly to capture various patterns of gene-gene interaction.
sis is such a simple combination of marginal regression and following penalized multiple regression  <cit>  that can become just an effective framework for the analysis of snp-snp interaction, which forms a typical ultrahigh-dimensional data in multiple regression model, regression coefficients can be used for evaluation of the importance of predictors, i.e.
coefficients closer to zero indicate less contribution to the regression formula.
statistical evaluation that infers regression coefficients to be zero is referred to as the variable selection  <cit> , which may lead to multiple regression model consisting of fewer predictors with improved estimation accuracy.
recent progress has justified theoretically and empirically the usefulness of variable selection in the large p small n situation  <cit> .
wu et al.
<cit> , hoggart et al.
<cit> , and ayers and cordell  <cit>  have been proposed applications of the modern variable selection to the selection of effective snps in single snp association study.
actually, sis has already been used in single snp association analysis in gwas data  <cit> .
recently, sis is theoretically proven to accept exponential grows of parameters as the increase of sample sizes in generalized linear models  <cit> .
this capability of exponentially large number of predictors relative to sample size encourages us the advanced application of the sis to the genome-wide gene-gene interaction analysis.
in its application, there is a difficulty in defining predictors for sis ranking step that can represent the interaction effect for each pair because the combination of two snps forms a  <dig> by  <dig> contingency table for which multiple patterns of interaction models are conceivable.
the ranking step is very important to effectively capture various patterns of interactions.
in this paper we elaborate the ranking strategy by proposing a promising three dummy coding methods and following variable selection procedure which is suitable for snp-snp interaction analysis.
we also implemented the proposed method in a very fast program, episis.
this is the first software program to employ the ultrahigh-dimensional variable selection method that can provide a statistically valid and high-speed exhaustive snp-snp interaction analysis even for standard gwas.
the application studies to simulated datasets show that the proposed method works successfully and accurately.
episis was applied to find some novel genetic components in the wtccc  data.
RESULTS
simulation experiments
we carried out simulation experiments to examine power and type  <dig> error rate for the proposed episis using the same dataset as those used in wan et al.
<cit>  for boost.
power
to investigate the power, we use the simulated data available from the boost website , which allows us a power comparison between episis, boost, and plink--fast-epistasis.
details of the datasets are summarized in additional file 1: text s <dig>  here we note that the power comparison between boost and plink--fast-epistasis has already been given by  <cit> .
the datasets were simulated from  <dig> different interaction models  based on the four epistatic interaction patterns  and three different mafs,  <dig> ,  <dig> , and  <dig> .
first three of the  <dig> scenarios are generated from model  <dig> with three mafs , next three are from model  <dig> , and so on.
each dataset contains  <dig> replicates with  <dig>  snps in 400/ <dig> or 800/ <dig> case–control individuals, among which the pair of first and last snps is set to the disease-susceptible combination and remaining  <dig> snps are non-risk factors.
power is calculated as the proportion that the interaction pair is detected in  <dig> replicates.
figures  <dig> and  <dig>  show the power results for three coding methods, likelihood cell-wise dummy coding , p-value cell-wise dummy coding , and p-value adaptive dummy coding , where the abscissa corresponds to the extended bayesian information criterion  tuning parameters γ between  <cit>  at an interval of  <dig> , i.e.
<dig> equally spaced points.
for comparison we applied the boost and plink –fast-epistasis for the same datasets.
the power results from these methods with significance levels of 5% and 30% are also shown in figures  <dig> and  <dig>  calculated by the output of boost’s p-value applying bonferroni correction on the basis of the number of interactions based on the number of snps, namely, we use the bonferroni correction with s/ <dig> hypotheses for s snps.
the bic  results in most powerful, then the power declines asγ goes to  <dig>  later, we discuss about the power carefully in terms of the balance with the type  <dig> error rate.
the  <dig> panels in figures  <dig> and  <dig> are arranged in the same order of the figure  <dig> of  <cit>  for ease of comparison.
because no versatile calibration method for γ has been proposed so far, we use the type  <dig> error rates resulted from simulated datasets as a guide for reasonable choice of γ.
type  <dig> error rate
we calculate the type  <dig> error rate as 1- .
we use the two scenarios given in  <cit> .
scenario 1: we generated  <dig>  replicates where no ld exists among snps using plink --simulate option, each having  <dig>  snps whose mafs are uniformly distributed within the interval  and including  <dig> case and  <dig> control individuals.
scenario 2: we generated  <dig> replicates where ld exists among snps using genomesimla  <cit>  on the basis of the marker information on the affymetrix 500Â k chip from human chromosome  <dig>  where each dataset contains  <dig>  snps in  <dig> case and  <dig> control samples.
the resulting type  <dig> error rates for episis are summarized in table  <dig> and those for plink –fast-epistasis and boost given in additional file 2: table s <dig> are close to the results in  <cit> .
in all simulations, bic  shows large type  <dig> error, and ebic at some γ >  <dig> leads to no type  <dig> error for all methods.
in summary, there exists a trade-off between power and type  <dig> error rate.
from the type  <dig> error rates in scenario  <dig>  corresponding to the null version of the power simulations, we can select  <dig>  for lcdc and pcdc, and  <dig>  for padc, at which the type  <dig> error rate is roughly comparable to the nominal error rates 10-30% adopted in figure  <dig> of  <cit> .
turning to the power plots at these ebic values, the padc shows highest power compared with the lcdc and pcdc, and outperforms boost except for cases  <dig>   <dig>  and  <dig> in 400/ <dig> samples, and in cases  <dig>   <dig>  and  <dig> for 800/ <dig> samples, which implies that the episis works effectively on the models that the boost shows low power.
on the other hand for the models on which the episis results in inferior performance to boost, the observation that the difference from boost in power is modest suggests that our proposed method can work effectively for several kinds of interaction patterns.
the optimal ebic values in scenario  <dig> were larger by  <dig>  than that in scenario  <dig> for three coding methods, lcdc, pcdc, and padc, which may come from the existence of ld or the fact that the number of snps before screening in scenario  <dig> is larger than that in scenario  <dig>
type  <dig> error rates in scenarios  <dig>  and  <dig>  with varying γ.
although it is impossible to complete simulation experiments in a genome-wide scale within realistic computational time, we adopt γ to be  <dig>  or  <dig>  for lcdc and pcdc, and  <dig>  or  <dig>  for padc, in the following wtccc data analysis, so that both power and type  <dig> error rate are of practical use as suggested by our simulation studies.
wtccc data analysis
we applied our program episis to real gwas datasets provided from wtccc , which included about  <dig>  common snp genotypes for each  <dig>  cases of seven human diseases, bipolar disorder , coronary artery disease , crohn’s disease , hypertension , rheumatoid arthritis , type  <dig> and type  <dig> diabetes , with  <dig>  controls  <cit> .
we first made high quality snp datasets through a standard quality-control filter, maf in control >  <dig> , hardy-weinberg exact  test p-value <  <dig> e- <dig>  study-wise missing data proportion >  <dig> , 1df trend test or 2df general test p-values <  <dig> e- <dig> between two control populations.
after eliminating snps without flag of “good-clustering” on signal summary information, for these diseases, we finally had datasets consisting of  <dig>  snp genotypes that were processed in our program.
for an exhaustive search of all possible two-pair interactions in addition to their main effects, on average, we needed about seven hours to complete one search on a standard linux machine having four gpu units .
we carried out three search strategies described in method section for each dataset using three different combinations of two ranking measures and two coding methods, i.e.
1) cdc with likelihood , 2) cdc with p-value , and 3) adc with p-value .
the details are described in method section.
we gathered the annotated information about each snp in our results from annotations released by affymetrix .
as summarized in table  <dig>  through these exhaustive searches, we finally obtained 1) seven interactions for bd, one for cad, one for cd,  <dig> for ht,  <dig> for ra, two for t1d and zero for t2d using lcdc , 2) two for bd, zero for cad, one for cd, two for ht, one for ra, one for t1d and zero for t2d using pcdc , 3) zero for bd, zero for cad, five for cd, zero for ht, one for ra, one for t1d and one for t2d using padc .
all of our results are listed in additional file 3: table s <dig>  among surviving variables through these searches, one main effect  for ra and two  for t1d were detected by lcdc as well as pcdc, which were found in the original report  <cit> .
in contrast, padc reported no main effect that was stronger than interactions.
we found a large number of interactions between snps located within a single gene locus  particularly in bd, ht and ra when using lcdc.
aepistasis; bdominance effect between snps in narrow  region or same locus; cmain effect; dtotal number of second-order interactions; enetwork group of interactions.
we detected the greatest number of interactions in ht by lcdc, bd and ht by pcdc, and cd by padc.
however, they included many redundant interactions that was shared a same snp as a partner.
for example, although padc reported five interactions in cd , those interactions shared same snps that can be assembled into two “network groups” implying potential higher-order interactions .
three methods in episis, lcdc, pcdc, and padc, yielded apparently quite different results in number of detected interactions, but by expressing as gene networks it turns out that they report similar gene-gene interactions in addition to some otherwise extras.
for example, in cd , some interesting interactions including a gene showing main effect were reported; that between atg16l <dig>  on 2q <dig>  and pdyn on 20p <dig>  =  <dig>  ) in addition to atg16l <dig> - ptger <dig>  - znf <dig>   = 2; ptger <dig> – znf300: oradc =  <dig> ) and so on.
lcdc appears to have a tendency to report a greater number of interactions than other two methods that utilize p-value ranking, and the difference was emphasized in dominance effects.
throughout our data analysis, the p-value ranking tends to be more conservative than the likelihood ranking.
on the other hand, cdcs appear to have a tendency to report dominance effects in addition to epistatic effects.
this phenomenon may be accounted for by the fact that the dominance effects are observed between linked loci, and so increases the sparsity of the  <dig> by  <dig> genotype contingency table because of the concentration to the diagonal cells .
a typical one of the examples is the second snp-snp pair found in the result for bd with the lcdc  given in additional file 3: table s <dig> online, where an apparent difference between cases and controls captured by the lcdc is seen in the aa-bb cell while other cells have very small number of observations due to the increased sparstiy.
as implied by the power simulations, cdc can work well for capturing disease-development models that suffice with a single cell, and therefore is applicable to the dominance effect mentioned above.
in cd, most of the interactions shared a snp in atg16l <dig>  which is an autophagy related gene showing the strongest main effect  originally reported in wtccc paper.
it is not always necessary that genes involving in interactions show strong main effect.
in another type of disease, bd, no strong main effect was reported, but episis detected a substantial number of interactions .
range of odds ratio  detected was between  <dig>  and  <dig> .
lcdc and pcdc reported a gene-gene interaction around mrpl <dig> .
CONCLUSIONS
we proposed an effective method for gene-gene interaction analysis using snps and developed a fast software program episis enabling genome-wide gene-gene interaction analysis by utilizing the cost-effective gpgpu technology.
this is the first method that successfully implements ultrahigh-dimensional variable selection approach for an exhaustive search for gene-gene interactions using realistic scale of snp-gwas data.
the method implemented employs a framework of sis  <cit> , which enables us to handle huge set of snp-snp interactions based on the modern large p small n regression methodology, rather than prohibitively conservative methods through conventional hypothesis testing.
simulation studies describe that our episis show successful performance.
among three ranking strategies proposed, padc showed greater performance than lcdc and pcdc in terms of power for most scenarios.
the exceptions exist in models  <dig> with maf =  <dig>   and model  <dig> with maf =  <dig>  , which may be well approximated by the disease-development model that suffices with a single cell of the  <dig> by  <dig> genotype interaction table due to the low mafs.
indeed, the cdc frequently captured the major homozogotes pair  in both two simulations.
intuitively, the cdc is workable for the interaction pattern influenced by a single or a few cells, while the adc is appealing to more complicated interaction patterns since the adc strengthens summary statistic of interaction by gathering cell-wise dummy predictors of cdc in an adaptive manner.
following the power simulation studies, we recommend using padc as a main method and cdc as a useful complement because the cdc may be of help in certain situations.
thus, we recommend using padc as a main method and cdcs as complements.
the power analysis through simulated datasets reveals that padc captures various interaction patterns, including the models for which boost resulted in low power, in practically high performance.
on the other hand, there exist some models where padc shows modestly lower power than boost but still maintains practicable performance.
our studies therefore show that the padc can have practically high power for several diseases models, which is a result from the use of analogous technique used in mdr by classifying snp-snp genotype pairs in two groups .
notably, since our proposed strategies generate  <dig> by  <dig> contingency table based on dummy coding, the sparse cell issue present in the boost is resolved.
one remaining issue in episis is that the control of type  <dig> error rates is somewhat rough, which is common in variable selection approaches, although it brings flexibility of the method.
the proposed choice of ebic parameter depends on the simulated datasets , which are less than genome-wide scale.
in addition, we observed modest increase of type  <dig> error rates under presence of linkage disequilibrium  compared with the simulation under no ld.
this observation is explained by a well-known fact that stringent correlation worsens performance of multiple regression.
nevertheless, additional file 5: figure s <dig> shows that lcdc and pcdc at γ =  <dig>  or padc at γ =  <dig>  in the application to seven real wtccc data can reduce the number of detected predictors as seen in the simulation studies, which implies that these ebic values for each screening method become a realistic solution.
here we note that the information about the number of detected predictors as varying γ may give us a further candidate gene sets.
indeed, the power plots imply that the smaller γ successfully increases the chance to include genuine interaction in the sis ranking.
we therefore encourage calibrating γ to be able to have the interaction pairs ranked high on the list, which could include susceptible loci to be examined in confirmatory study.
such a flexible usage is one of the advantages in using machine learning-based methods.
by applying episis to wtccc datasets, we obtained a large number of interactions potentially conferring susceptibility to them.
the distributions of the interactions detected episis were disease-specific but not software- or method-specific, implying that these results were likely derived due to the genuine distributions of these interaction but not pseudo-negative and/or positive due to the algorithm in episis.
although a confirmatory study using independent sample set is needed to eliminate bias and confounding specific to original data, some of these interactions assembled in an interpretable network graph appear to be plausible from their functional points of view.
for example, a statistical interaction among atg16l1/ptger4/znf <dig> in cd could imply the involvement of a synergistic combination among autophagy and other possible mechanisms such as immune-inflammation in the etiology of inflammatory bowel disease, which have been repeatedly suggested by genetic studies about their main effects  <cit>  as well as mechanistic studies  <cit> .
in contrast, one network of interactions found in bd centered mrpl <dig> connected to many other genes of which some are possibly involved in mitochondrial function.
this finding seems to support the previous hypothesis that mitochondrial functions underlies the pathophysiology of mental illnesses such as bipolar disorder and schizophrenia  <cit> .
episis also reported a large number of interactions between neighboring snps, which are expected to have remarkable dominance effects through potential haplotype-blocks, which might confer disease susceptibility as haplotype blocks consisted of snps in a single locus.
for illustrative purpose, we examined the ld pattern of the dominance effect detected by episis using two genes, ppm1a and ulk <dig>  as an example, which were found in our analysis for ht.
as stated above, the episis attempts to detect the difference between case- and control-distributions.
to see this in detail, we provide in figure  <dig> the correlation coefficients based on haplotype and genotype frequencies, which can represent a distributional characteristic.
the genotype-based correlation coefficient is proposed by wellek and ziegler  <cit> .
to compare the difference between cases and controls, wellek and ziegler’s statistic is more appropriate than the haplotype-based correlation coefficient, since the estimation of haplotype frequencies requires the condition in which hardy-weinberg equilibrium holds, which is not expected in case population.
since the purpose here is to compare case- and control-distributions, consideration of the sign is needed.
figure  <dig> describes the correlation coefficient matrix for cases and controls for ppm1a  and ulk <dig> .
haplotype-based and genotype-based correlations are given in the upper- and lower-triangle parts, respectively.
snp-snp pairs that episis detected are emphasized by red diamond symbol.
it can be seen that all pairs are included in pairs highlighted in the figure.
here we compare our results in wtccc data analysis with those using boost and plink –fast-epistasis respectively given by wan et al.
<cit>  and cordell  <cit> .
first, since exhaustive search using plink –fast-epistasis is virtually infeasible due to computational burden, cordell  <cit>  conducts semi-exhaustive search for cd based on snps that passed the single locus p-value threshold of  <dig> .
she concluded that the snp-snp pair that showed highly significant fast-epistasis p-value is false positive.
on the other hand, episis detected a few interactions for cd in lcdc, pcdc, and padc as shown in additional file 4: table s <dig>  the analysis using episis contains no interactions located in top  <dig> using plink given in additional file 2: table s <dig> of  <cit> .
this difference comes from the semi-exhaustive search used in  <cit>  as well as the fact that the summary statistic in fast-epistasis is based on allelic correlation whereas that in episis is based on genotypic data.
second, wan et al.
<cit>  report many gene-gene interactions detected by the boost for t1d in the mhc  region, although they mention that the interactions found for other six diseases are nontrivial except for one snp pair in cd .
they state that the interactions obtained in t1d analysis after excluding significant loci under single-locus scan include interesting interaction patterns between the mhc class i and class ii.
although our episis finds some interactions at two snps in the mhc class ii, rs <dig> and rs <dig>  these two snps show strong main effect and are interacted with snps on other chromosome, i.e.
not in the mhc class i. this discrepancy from the boost’s result could be partly explained by the fact that we analyzed entire snps without imposing single-locus significance threshold.
although it could be possible to exclude snps showing strong main effect, we prefer including them since they contribute forming the affected and unaffected populations and could be ancillary to interaction analysis through regression model.
the episis’s result on t1d tells that the snps in the mhc class ii region can play a major role in partitioning affected and unaffected individuals even when considering second-order interactions.
it is also noteworthy that the interaction detected by episis include various snp-snp pairs having cells with fewer observations , implying that the episis overcomes the sparse cell issue.
since boost, plink –fast-epistasis, and episis have advantage and disadvantage depending on disease patterns that underlie as seen in power simulation, we recommend using them in a mutually complementary manner.
finally, we summarize our future works in what follows.
first, since we use the multiple regression model, inclusion of covariates is possible to adjust the influence of several confounding factors, such as sex, age, or population stratification.
second, the generalized linear model used in sis allows extension to more general phenotypes other than disease status.
third, although the current version does not allow the presence of missing genotype in the future version, we would overcome this issue.
fourth, an iterative version of sis  has been proposed  <cit> , which could improve the detection ability.
we have already implemented the isis in episis, which is useful in finding additional interactions confounded by interactions that have already detected.
it poses further computational cost due to increased dimensionality of parameters, where the cost depends on the number of detected factors in sis and subsequent penalized regression.
an issue arises when detected dummy predictors after sis are highly-correlated, i.e.
multi-collinearlity, causing failure in convergence of estimating regression coefficients in sis ranking step, which incurs considerable increase of computational time while some of results lose their reliability.
unfortunately, we provide no reasonable solution to this issue and it is still under consideration.
isis could work if no multi-collinearlity exists, although we have no knowledge about relationship between type  <dig> error rates and ebic tuning parameters, with extra predictors provided from the previous sis step, in penalized regression step.
consequently, at the moment, we have to pay adequate attention to the multi-collinearity within the dummy predictors as well as to the selection of ebic tuning parameter in using isis procedure, and we will address such issues in our future work.
although the isis where the above issues are resolved is ultimately desirable, the current episis will work reasonably as in simple situations used in the simulation studies.
our program episis described in this manuscript will be freely available from the authors’ webpage soon.

 single snp association study is a popular method to detect genes that are susceptible to human diseases. although candidate gene approach that uses prior knowledge about its function is an efficient procedure, it could overlook genes whose function is unknown or ambiguous. gwas using whole genome snps is thus an attractive solution to this issue. many procedures proposed for gwas so far are based on marginal association between each snp and phenotype. however, it has been pointed out that most susceptible-snps identified often show low or moderate effect size, and hence may explain only a few percentage of the genetic variance. the fact that, for several complex diseases, the recurrence risk ratio decreases quicker than 1/ as the relatedness decreases implies the involvement of non-additive interactions in their etiology. in addition to substantial contributions of rarer variants, gene-gene and gene-environment interaction would be one of strong candidates that can explain the missing heritability as well. without investigating such interactions, we therefore may miss genuine disease-susceptible loci. an effective and accurate method to search gene-gene interactions will utilize immediately original snp-gwas data to decipher such missing genetic components of complex human diseases. in this paper, we tackle a development of powerful method for the genome-wide gene-gene interaction analysis using snps. marchinni et al. suggest that the use of arbitrary single locus-disease association model under which some interaction effect is present might prevent from finding susceptible interaction effect. therefore exhaustive search for interactions are needed rather than primary screening by marginal effect. although individual hypothesis tests for all snp-snp pairs as frequently used in snp-gwas might be the simplest approach, it is however faced with multiple testing issues caused by -, snps. specifically, there is a difficulty in setting genome-wide significance level using, e.g., bonferroni correction, fprp, or bayes factor, leading to prohibitively conservative result because they fail to successfully incorporate the correlation structure between each hypothesis. for instance, the total number of hypothesises for second-order gene-gene interaction is about – in standard gwas data. then bonferroni-corrected significance level must be considerably small by multiplying the nominal significance level such as by the correction factor less than 10− no efficient and universal multiple testing method to deal with such the huge set of hypotheses having complicated correlation structure is proposed so far. although multiple logistic regression may be the method that can incorporate the correlation structure between hypothesises by putting them as predictor variables, the so-called large p small n situation, precludes simultaneous inclusion of all snp-snp interactions in the logistic regression model, namely, no unique solution is determined. there are several software programs for gene-gene interaction analysis including multifactor dimensionality reduction, fast-epistasis option in plink, tuning relief, random jungle, beam, and boost. cordell carried out an extensive comparative study for these methods. she reported that mdr and beam have computational difficulty in analyzing current scale in gwas and that random jungle is applicable to small scale dataset but not for whole-genome dataset. plink-fast-epistasis and boost are designed for exhaustive search. since both methods are based on hypothesis testing principles, as stated above, exact determination of significance level for whole-genome data is difficult due to the numerous multiple hypotheses in addition to the complicated correlation structure between them. hypothesis testing methods could be too conservative and lead to false negative result, because some complicated correlation structure would commonly appear when considering interaction pairs. the boost has been shown to outperform the plink--fast-epistasis in most of the interaction models considered in. the boost primarily performs screening of possible hypotheses using kirkwood superposition approximation, then applies standard multiple testing on the basis of the initial number of snp-snp pairs to the likelihood ratio statistic under the saturated logistic regression model having the degrees of freedom of however, the test statistic may fail to follow the chi-squared distribution of df under the presence of sparse cells, which causes too small type error rate. because the sparse cell issue is unavoidable in genotypic interaction data, some solutions in development of summary statistic for gene-gene interaction are strongly necessitated. our method is free from the sparse cell issue owing to the proposed dummy coding strategy, which generates by contingency table. in this paper we develop a new method based on contemporary ultrahigh-dimensional variable selection approach, designated the sure independence screening rather than hypothesis testing. the principle of sis is based on the same as that of machine learning methods exploring the optimal model that can explain sufficiently the data in the most parsimonious way; its methodological feature to avoid over- or under-fitting enables us flexibly to capture various patterns of gene-gene interaction. sis is such a simple combination of marginal regression and following penalized multiple regression that can become just an effective framework for the analysis of snp-snp interaction, which forms a typical ultrahigh-dimensional data in multiple regression model, regression coefficients can be used for evaluation of the importance of predictors, i.e. coefficients closer to zero indicate less contribution to the regression formula. statistical evaluation that infers regression coefficients to be zero is referred to as the variable selection, which may lead to multiple regression model consisting of fewer predictors with improved estimation accuracy. recent progress has justified theoretically and empirically the usefulness of variable selection in the large p small n situation. wu et al., hoggart et al., and ayers and cordell have been proposed applications of the modern variable selection to the selection of effective snps in single snp association study. actually, sis has already been used in single snp association analysis in gwas data. recently, sis is theoretically proven to accept exponential grows of parameters as the increase of sample sizes in generalized linear models. this capability of exponentially large number of predictors relative to sample size encourages us the advanced application of the sis to the genome-wide gene-gene interaction analysis. in its application, there is a difficulty in defining predictors for sis ranking step that can represent the interaction effect for each pair because the combination of two snps forms a by contingency table for which multiple patterns of interaction models are conceivable. the ranking step is very important to effectively capture various patterns of interactions. in this paper we elaborate the ranking strategy by proposing a promising three dummy coding methods and following variable selection procedure which is suitable for snp-snp interaction analysis. we also implemented the proposed method in a very fast program, episis. this is the first software program to employ the ultrahigh-dimensional variable selection method that can provide a statistically valid and high-speed exhaustive snp-snp interaction analysis even for standard gwas. the application studies to simulated datasets show that the proposed method works successfully and accurately. episis was applied to find some novel genetic components in the wtccc data.  simulation experiments we carried out simulation experiments to examine power and type error rate for the proposed episis using the same dataset as those used in wan et al. for boost. power to investigate the power, we use the simulated data available from the boost website, which allows us a power comparison between episis, boost, and plink--fast-epistasis. details of the datasets are summarized in additional file 1: text s here we note that the power comparison between boost and plink--fast-epistasis has already been given by. the datasets were simulated from different interaction models based on the four epistatic interaction patterns and three different mafs,,, and. first three of the scenarios are generated from model with three mafs, next three are from model, and so on. each dataset contains replicates with snps in 400/ or 800/ case–control individuals, among which the pair of first and last snps is set to the disease-susceptible combination and remaining snps are non-risk factors. power is calculated as the proportion that the interaction pair is detected in replicates. figures and show the power results for three coding methods, likelihood cell-wise dummy coding, p-value cell-wise dummy coding, and p-value adaptive dummy coding, where the abscissa corresponds to the extended bayesian information criterion tuning parameters γ between at an interval of, i.e. equally spaced points. for comparison we applied the boost and plink –fast-epistasis for the same datasets. the power results from these methods with significance levels of 5% and 30% are also shown in figures and calculated by the output of boost’s p-value applying bonferroni correction on the basis of the number of interactions based on the number of snps, namely, we use the bonferroni correction with s/ hypotheses for s snps. the bic results in most powerful, then the power declines asγ goes to later, we discuss about the power carefully in terms of the balance with the type error rate. the panels in figures and are arranged in the same order of the figure of for ease of comparison. because no versatile calibration method for γ has been proposed so far, we use the type error rates resulted from simulated datasets as a guide for reasonable choice of γ. type error rate we calculate the type error rate as 1-. we use the two scenarios given in. scenario 1: we generated replicates where no ld exists among snps using plink --simulate option, each having snps whose mafs are uniformly distributed within the interval and including case and control individuals. scenario 2: we generated replicates where ld exists among snps using genomesimla on the basis of the marker information on the affymetrix 500Â k chip from human chromosome where each dataset contains snps in case and control samples. the resulting type error rates for episis are summarized in table and those for plink –fast-epistasis and boost given in additional file 2: table s are close to the results in. in all simulations, bic shows large type error, and ebic at some γ > leads to no type error for all methods. in summary, there exists a trade-off between power and type error rate. from the type error rates in scenario corresponding to the null version of the power simulations, we can select for lcdc and pcdc, and for padc, at which the type error rate is roughly comparable to the nominal error rates 10-30% adopted in figure of. turning to the power plots at these ebic values, the padc shows highest power compared with the lcdc and pcdc, and outperforms boost except for cases and in 400/ samples, and in cases and for 800/ samples, which implies that the episis works effectively on the models that the boost shows low power. on the other hand for the models on which the episis results in inferior performance to boost, the observation that the difference from boost in power is modest suggests that our proposed method can work effectively for several kinds of interaction patterns. the optimal ebic values in scenario were larger by than that in scenario for three coding methods, lcdc, pcdc, and padc, which may come from the existence of ld or the fact that the number of snps before screening in scenario is larger than that in scenario type error rates in scenarios and with varying γ. although it is impossible to complete simulation experiments in a genome-wide scale within realistic computational time, we adopt γ to be or for lcdc and pcdc, and or for padc, in the following wtccc data analysis, so that both power and type error rate are of practical use as suggested by our simulation studies. wtccc data analysis we applied our program episis to real gwas datasets provided from wtccc, which included about common snp genotypes for each cases of seven human diseases, bipolar disorder, coronary artery disease, crohn’s disease, hypertension, rheumatoid arthritis, type and type diabetes, with controls. we first made high quality snp datasets through a standard quality-control filter, maf in control >, hardy-weinberg exact test p-value < e- study-wise missing data proportion >, 1df trend test or 2df general test p-values < e- between two control populations. after eliminating snps without flag of “good-clustering” on signal summary information, for these diseases, we finally had datasets consisting of snp genotypes that were processed in our program. for an exhaustive search of all possible two-pair interactions in addition to their main effects, on average, we needed about seven hours to complete one search on a standard linux machine having four gpu units. we carried out three search strategies described in method section for each dataset using three different combinations of two ranking measures and two coding methods, i.e. 1) cdc with likelihood, 2) cdc with p-value, and 3) adc with p-value. the details are described in method section. we gathered the annotated information about each snp in our results from annotations released by affymetrix. as summarized in table through these exhaustive searches, we finally obtained 1) seven interactions for bd, one for cad, one for cd, for ht, for ra, two for t1d and zero for t2d using lcdc, 2) two for bd, zero for cad, one for cd, two for ht, one for ra, one for t1d and zero for t2d using pcdc, 3) zero for bd, zero for cad, five for cd, zero for ht, one for ra, one for t1d and one for t2d using padc. all of our results are listed in additional file 3: table s among surviving variables through these searches, one main effect for ra and two for t1d were detected by lcdc as well as pcdc, which were found in the original report. in contrast, padc reported no main effect that was stronger than interactions. we found a large number of interactions between snps located within a single gene locus particularly in bd, ht and ra when using lcdc. aepistasis; bdominance effect between snps in narrow region or same locus; cmain effect; dtotal number of second-order interactions; enetwork group of interactions. we detected the greatest number of interactions in ht by lcdc, bd and ht by pcdc, and cd by padc. however, they included many redundant interactions that was shared a same snp as a partner. for example, although padc reported five interactions in cd, those interactions shared same snps that can be assembled into two “network groups” implying potential higher-order interactions. three methods in episis, lcdc, pcdc, and padc, yielded apparently quite different results in number of detected interactions, but by expressing as gene networks it turns out that they report similar gene-gene interactions in addition to some otherwise extras. for example, in cd, some interesting interactions including a gene showing main effect were reported; that between atg16l on 2q and pdyn on 20p = ) in addition to atg16l - ptger - znf = 2; ptger – znf300: oradc = ) and so on. lcdc appears to have a tendency to report a greater number of interactions than other two methods that utilize p-value ranking, and the difference was emphasized in dominance effects. throughout our data analysis, the p-value ranking tends to be more conservative than the likelihood ranking. on the other hand, cdcs appear to have a tendency to report dominance effects in addition to epistatic effects. this phenomenon may be accounted for by the fact that the dominance effects are observed between linked loci, and so increases the sparsity of the by genotype contingency table because of the concentration to the diagonal cells. a typical one of the examples is the second snp-snp pair found in the result for bd with the lcdc given in additional file 3: table s online, where an apparent difference between cases and controls captured by the lcdc is seen in the aa-bb cell while other cells have very small number of observations due to the increased sparstiy. as implied by the power simulations, cdc can work well for capturing disease-development models that suffice with a single cell, and therefore is applicable to the dominance effect mentioned above. in cd, most of the interactions shared a snp in atg16l which is an autophagy related gene showing the strongest main effect originally reported in wtccc paper. it is not always necessary that genes involving in interactions show strong main effect. in another type of disease, bd, no strong main effect was reported, but episis detected a substantial number of interactions. range of odds ratio detected was between and. lcdc and pcdc reported a gene-gene interaction around mrpl.  we proposed an effective method for gene-gene interaction analysis using snps and developed a fast software program episis enabling genome-wide gene-gene interaction analysis by utilizing the cost-effective gpgpu technology. this is the first method that successfully implements ultrahigh-dimensional variable selection approach for an exhaustive search for gene-gene interactions using realistic scale of snp-gwas data. the method implemented employs a framework of sis, which enables us to handle huge set of snp-snp interactions based on the modern large p small n regression methodology, rather than prohibitively conservative methods through conventional hypothesis testing. simulation studies describe that our episis show successful performance. among three ranking strategies proposed, padc showed greater performance than lcdc and pcdc in terms of power for most scenarios. the exceptions exist in models with maf = and model with maf =, which may be well approximated by the disease-development model that suffices with a single cell of the by genotype interaction table due to the low mafs. indeed, the cdc frequently captured the major homozogotes pair in both two simulations. intuitively, the cdc is workable for the interaction pattern influenced by a single or a few cells, while the adc is appealing to more complicated interaction patterns since the adc strengthens summary statistic of interaction by gathering cell-wise dummy predictors of cdc in an adaptive manner. following the power simulation studies, we recommend using padc as a main method and cdc as a useful complement because the cdc may be of help in certain situations. thus, we recommend using padc as a main method and cdcs as complements. the power analysis through simulated datasets reveals that padc captures various interaction patterns, including the models for which boost resulted in low power, in practically high performance. on the other hand, there exist some models where padc shows modestly lower power than boost but still maintains practicable performance. our studies therefore show that the padc can have practically high power for several diseases models, which is a result from the use of analogous technique used in mdr by classifying snp-snp genotype pairs in two groups. notably, since our proposed strategies generate by contingency table based on dummy coding, the sparse cell issue present in the boost is resolved. one remaining issue in episis is that the control of type error rates is somewhat rough, which is common in variable selection approaches, although it brings flexibility of the method. the proposed choice of ebic parameter depends on the simulated datasets, which are less than genome-wide scale. in addition, we observed modest increase of type error rates under presence of linkage disequilibrium compared with the simulation under no ld. this observation is explained by a well-known fact that stringent correlation worsens performance of multiple regression. nevertheless, additional file 5: figure s shows that lcdc and pcdc at γ = or padc at γ = in the application to seven real wtccc data can reduce the number of detected predictors as seen in the simulation studies, which implies that these ebic values for each screening method become a realistic solution. here we note that the information about the number of detected predictors as varying γ may give us a further candidate gene sets. indeed, the power plots imply that the smaller γ successfully increases the chance to include genuine interaction in the sis ranking. we therefore encourage calibrating γ to be able to have the interaction pairs ranked high on the list, which could include susceptible loci to be examined in confirmatory study. such a flexible usage is one of the advantages in using machine learning-based methods. by applying episis to wtccc datasets, we obtained a large number of interactions potentially conferring susceptibility to them. the distributions of the interactions detected episis were disease-specific but not software- or method-specific, implying that these results were likely derived due to the genuine distributions of these interaction but not pseudo-negative and/or positive due to the algorithm in episis. although a confirmatory study using independent sample set is needed to eliminate bias and confounding specific to original data, some of these interactions assembled in an interpretable network graph appear to be plausible from their functional points of view. for example, a statistical interaction among atg16l1/ptger4/znf in cd could imply the involvement of a synergistic combination among autophagy and other possible mechanisms such as immune-inflammation in the etiology of inflammatory bowel disease, which have been repeatedly suggested by genetic studies about their main effects as well as mechanistic studies. in contrast, one network of interactions found in bd centered mrpl connected to many other genes of which some are possibly involved in mitochondrial function. this finding seems to support the previous hypothesis that mitochondrial functions underlies the pathophysiology of mental illnesses such as bipolar disorder and schizophrenia. episis also reported a large number of interactions between neighboring snps, which are expected to have remarkable dominance effects through potential haplotype-blocks, which might confer disease susceptibility as haplotype blocks consisted of snps in a single locus. for illustrative purpose, we examined the ld pattern of the dominance effect detected by episis using two genes, ppm1a and ulk as an example, which were found in our analysis for ht. as stated above, the episis attempts to detect the difference between case- and control-distributions. to see this in detail, we provide in figure the correlation coefficients based on haplotype and genotype frequencies, which can represent a distributional characteristic. the genotype-based correlation coefficient is proposed by wellek and ziegler. to compare the difference between cases and controls, wellek and ziegler’s statistic is more appropriate than the haplotype-based correlation coefficient, since the estimation of haplotype frequencies requires the condition in which hardy-weinberg equilibrium holds, which is not expected in case population. since the purpose here is to compare case- and control-distributions, consideration of the sign is needed. figure describes the correlation coefficient matrix for cases and controls for ppm1a and ulk. haplotype-based and genotype-based correlations are given in the upper- and lower-triangle parts, respectively. snp-snp pairs that episis detected are emphasized by red diamond symbol. it can be seen that all pairs are included in pairs highlighted in the figure. here we compare our results in wtccc data analysis with those using boost and plink –fast-epistasis respectively given by wan et al. and cordell. first, since exhaustive search using plink –fast-epistasis is virtually infeasible due to computational burden, cordell conducts semi-exhaustive search for cd based on snps that passed the single locus p-value threshold of. she concluded that the snp-snp pair that showed highly significant fast-epistasis p-value is false positive. on the other hand, episis detected a few interactions for cd in lcdc, pcdc, and padc as shown in additional file 4: table s the analysis using episis contains no interactions located in top using plink given in additional file 2: table s of. this difference comes from the semi-exhaustive search used in as well as the fact that the summary statistic in fast-epistasis is based on allelic correlation whereas that in episis is based on genotypic data. second, wan et al. report many gene-gene interactions detected by the boost for t1d in the mhc region, although they mention that the interactions found for other six diseases are nontrivial except for one snp pair in cd. they state that the interactions obtained in t1d analysis after excluding significant loci under single-locus scan include interesting interaction patterns between the mhc class i and class ii. although our episis finds some interactions at two snps in the mhc class ii, rs and rs these two snps show strong main effect and are interacted with snps on other chromosome, i.e. not in the mhc class i. this discrepancy from the boost’s result could be partly explained by the fact that we analyzed entire snps without imposing single-locus significance threshold. although it could be possible to exclude snps showing strong main effect, we prefer including them since they contribute forming the affected and unaffected populations and could be ancillary to interaction analysis through regression model. the episis’s result on t1d tells that the snps in the mhc class ii region can play a major role in partitioning affected and unaffected individuals even when considering second-order interactions. it is also noteworthy that the interaction detected by episis include various snp-snp pairs having cells with fewer observations, implying that the episis overcomes the sparse cell issue. since boost, plink –fast-epistasis, and episis have advantage and disadvantage depending on disease patterns that underlie as seen in power simulation, we recommend using them in a mutually complementary manner. finally, we summarize our future works in what follows. first, since we use the multiple regression model, inclusion of covariates is possible to adjust the influence of several confounding factors, such as sex, age, or population stratification. second, the generalized linear model used in sis allows extension to more general phenotypes other than disease status. third, although the current version does not allow the presence of missing genotype in the future version, we would overcome this issue. fourth, an iterative version of sis has been proposed, which could improve the detection ability. we have already implemented the isis in episis, which is useful in finding additional interactions confounded by interactions that have already detected. it poses further computational cost due to increased dimensionality of parameters, where the cost depends on the number of detected factors in sis and subsequent penalized regression. an issue arises when detected dummy predictors after sis are highly-correlated, i.e. multi-collinearlity, causing failure in convergence of estimating regression coefficients in sis ranking step, which incurs considerable increase of computational time while some of results lose their reliability. unfortunately, we provide no reasonable solution to this issue and it is still under consideration. isis could work if no multi-collinearlity exists, although we have no knowledge about relationship between type error rates and ebic tuning parameters, with extra predictors provided from the previous sis step, in penalized regression step. consequently, at the moment, we have to pay adequate attention to the multi-collinearity within the dummy predictors as well as to the selection of ebic tuning parameter in using isis procedure, and we will address such issues in our future work. although the isis where the above issues are resolved is ultimately desirable, the current episis will work reasonably as in simple situations used in the simulation studies. our program episis described in this manuscript will be freely available from the authors’ webpage soon. 
episis can complete exhaustive search for snp-snp interactions in standard gwas dataset within several hours. genome-wide gene-gene interaction analysis using single nucleotide polymorphisms is an attractive way for identification of genetic components that confers susceptibility of human complex diseases. a large number of snp-snp pairs than sample size, so-called the large p small n problem, precludes simultaneous analysis using multiple regression. individual hypothesis testing for snp-snp pairs as in common genome-wide association study however involves difficulty in setting overall p-value due to complicated correlation structure, namely, the multiple testing problem that causes unacceptable false negative results. we adopt an up-to-date method for ultrahigh-dimensional variable selection termed the sure independence screening for appropriate handling of numerous number of snp-snp interactions by including them as predictor variables in logistic regression. we also implemented the procedures in a software program, episis, using the cost-effective gpgpu technology. the proposed method works successfully in simulation experiments and in application to real wtccc data. we propose ranking strategy using promising dummy coding methods and following variable selection procedure in the sis method suitably modified for gene-gene interaction analysis. based on the machine-learning principle, the proposed method gives powerful and flexible genome-wide search for various patterns of gene-gene interaction.
 genome-wide gene-gene interaction analysis using single nucleotide polymorphisms is an attractive way for identification of genetic components that confers susceptibility of human complex diseases. individual hypothesis testing for snp-snp pairs as in common genome-wide association study however involves difficulty in setting overall p-value due to complicated correlation structure, namely, the multiple testing problem that causes unacceptable false negative results. a large number of snp-snp pairs than sample size, so-called the large p small n problem, precludes simultaneous analysis using multiple regression. the method that overcomes above issues is thus needed.  we adopt an up-to-date method for ultrahigh-dimensional variable selection termed the sure independence screening for appropriate handling of numerous number of snp-snp interactions by including them as predictor variables in logistic regression. we propose ranking strategy using promising dummy coding methods and following variable selection procedure in the sis method suitably modified for gene-gene interaction analysis. we also implemented the procedures in a software program, episis, using the cost-effective gpgpu technology. episis can complete exhaustive search for snp-snp interactions in standard gwas dataset within several hours. the proposed method works successfully in simulation experiments and in application to real wtccc data.  based on the machine-learning principle, the proposed method gives powerful and flexible genome-wide search for various patterns of gene-gene interaction. 