A Novel Statistic for Genome-Wide Interaction Analysis

Although great progress in genome-wide association studies (GWAS) has been made, the significant SNP associations identified by GWAS account for only a few percent of the genetic variance, leading many to question where and how we can find the missing heritability. There is increasing interest in genome-wide interaction analysis as a possible source of finding heritability unexplained by current GWAS. However, the existing statistics for testing interaction have low power for genome-wide interaction analysis. To meet challenges raised by genome-wide interactional analysis, we have developed a novel statistic for testing interaction between two loci (either linked or unlinked). The null distribution and the type I error rates of the new statistic for testing interaction are validated using simulations. Extensive power studies show that the developed statistic has much higher power to detect interaction than classical logistic regression. The results identified 44 and 211 pairs of SNPs showing significant evidence of interactions with FDR<0.001 and 0.001<FDR<0.003, respectively, which were seen in two independent studies of psoriasis. These included five interacting pairs of SNPs in genes LST1/NCR3, CXCR5/BCL9L, and GLS2, some of which were located in the target sites of miR-324-3p, miR-433, and miR-382, as well as 15 pairs of interacting SNPs that had nonsynonymous substitutions. Our results demonstrated that genome-wide interaction analysis is a valuable tool for finding remaining missing heritability unexplained by the current GWAS, and the developed novel statistic is able to search significant interaction between SNPs across the genome. Real data analysis showed that the results of genome-wide interaction analysis can be replicated in two independent studies.


Introduction
In the past three years, about 400 genome-wide association studies (GWAS) that focused largely on individually testing the associations of single SNP with diseases have been conducted [1]. These studies have identified more than 531 SNPs associated with different traits or diseases [2] and have provided substantial information for understanding disease mechanisms. Despite the progress that has been made, the significant SNP associations identified by GWAS account for only a few percent of the genetic variance which begs the question where and how the missing heritability can be identified [3,4]. Possible explanations include [1,4]: (1) The previous GWAS are mainly based on the common disease, common variant hypothesis. However, in addition to single nucleotide polymorphisms (SNPs) with a minor allele frequency (MAF) greater than 1%, there are other classes of human genetic variation including: (a) rare variants that are defined as mutations with a MAF of less than 1% and (b) structural variants including copy number variants (CNVs) and copy neutral variation such as inversions and translocations. Common diseases can also be caused by multiple rare mutations, each with a low marginal genetic effect. A more realistic model is that the entire spectrum of genetic variants ranging from rare to common contributes to disease susceptibility. (2) Most of current GWAS have focused on SNP analysis in which each variant is tested for association individually. However, common disease often arises from the combined effect of multiple loci within a gene or interaction of multiple genes within a pathway. If we only consider the most significant SNPs, the genetic variants that jointly have significant impact on risk, but individually make only a small contribution, will be missed. Another way to discover the missing heritability of complex diseases is to investigate gene-gene and gene-environment interaction. Disease development is a dynamic process of genegene and gene-environment interactions within a complex biological system which is organized into interacting networks [5]. Modern complexity theory assumes that the complexity is attributed to the interactions among the components of the system, therefore, interaction has been considered as a sensible measure of complexity of the biological systems. The more interactions between the components there are, the more complex the system is. The disease may be caused by joint action of multiple loci. Motivation for studying statistical interaction is to provide increased power for detecting joint acting effects of interacting loci than testing for only marginal association of each of the loci individually. Screening for only main effects might miss the vast majority of the genetic variants that interact with each other and with environment to cause diseases [6]. We argue that the interactions hold a key for dissecting the genetic structure of complex diseases and elucidating the biological and biochemical pathway underlying the diseases [7,8]. Ignoring gene-gene and gene-environment interactions will likely obscure the detection of genetic effects and may lead to inconsistent association results across studies [9,10].
GWAS in which several hundred thousands or even a millions of SNPs are typed in thousands of individuals provide unprecedented opportunities for systematic exploration of the universe of variants and interactions in the entire genome and also raise several serious challenges for genome-wide interaction analysis. The first challenge comes from the problems imposed by multiple testing. Even for investigating pair-wise interaction, the total number of tests for interaction between all possible SNPs across the genome will be extremely large. Bonferroni-corrected P-values for ensuring genome-wide significance level of 0.05 will be too small to reach. The second challenge is the need for computationally simple statistics for testing interactions. The simplest way to search for interactions between two loci is to test all possible two-locus interactions. This exhaustive search demands large computations. Therefore, the computational time of each twolocus interaction test should be short. The third challenge is the power of the statistics for testing interaction. To ensure the genome-wide significance, the statistics should have high power to detect interaction. Developing simple and efficient analytic methods for evaluation of the gene-gene interactions is critical to the success of genome-wide gene-gene interaction analysis. Finally, the fourth challenge is replication of the finding of such interactions in independent studies. This report will attempt to meet these challenges, at least in part. To achieve this, we first should define a good measure of gene-gene interaction. Despite current enthusiasm for investigation of genegene interactions, published results that document these interactions in humans are limited and the essential issue of how to define and detect gene-gene interactions remains unresolved. Over the last three decades, epidemiologists have debated intensely about how to define and measure interaction in epidemiologic studies [7,8,[11][12][13][14][15]; The concept of gene-gene interactions is often used, but rarely specified with precision [16]. In general, statistical gene-gene interaction is defined as departure from additive or multiplicative joint effects of the genetic risk factors [17]. It is increasingly recognized that statistical interactions are scale dependent [18]. In other words, how to define the effects of a risk factor and how to measure departure from the independence of effects will greatly affect assessment of gene-gene interaction. The most popular scale upon which risk factors are measured in case-control studies is oddsratio. The traditional odds-ratio is defined in terms of genotypes at two loci. Similar to two-locus association analysis where only genotype information at two loci is used, odds-ratio defined by genotypes for testing interaction will not employ allelic association information. However, it is known that interaction between two loci will generate allelic associations in some circumstances [19]. Since they do not use allelic association information between two loci, the statistical methods based on the odds-ratio that is defined in terms of genotypes will have less power to detect interaction. To overcome this limitation, we will define odds-ratio in terms of a pseudohaplotype (which is defined as two alleles located on the same paternal or maternal chromosomes) for measuring interaction, and then we will investigate its properties and develop a statistic based on pseudohaplotype defined odds-ratio for testing interaction between two loci (either linked or unlinked).
To demonstrate that the pseudohaplotype odds-ratio interaction measure-based statistic for detection of interaction between two loci will not cause false positive problems, we then investigate type I error rates. To reveal the merit and limitation of the pseudohaplotype odds-ratio interaction measure-based statistic for detection of interaction, we will compare its power for detecting interaction with the traditional logistic regression and ''fast-epistasis'' in PLINK [20].
Although nearly 400 GWAS have been documented, few genome-wide interaction analyses have been performed and few findings of significant interaction reported [8,21,22]. Emily et al [23] tested about 3,107,904-3,850,339 pairs of SNPs located in genes with potential protein-protein interaction and reported four significant cases of interactions, one in each of Crohn's Disease, bipolar disorder, hypertension and rheumatoid arthritis in the WTCCC dataset, but these have not been replicated. To further evaluate the performance of our new statistic and test the feasibility of genome-wide interaction analysis, the presented statistic was applied to interaction analysis of two independent GWAS datasets of psoriasis where 1,266,378,301 pairs of SNPs from 50,327 SNPs in the first dataset and 1,243,782,750 pairs of SNPs from 49,876 SNPs in the second dataset were tested for interactions. These SNPs in the datasets were selected from 501 pathways assembled from KEGG [24] and Biocarta (http://www. biocarta.com) pathway databases. A program for using the developed statistic to test interaction which was implemented by C++ can be downloaded from our website http://www.sph.uth. tmc.edu/hgc/faculty/xiong/index.htm.

Methods
A case-control study design for detection of interaction between two loci (SNPs) where two loci can be either linked or unlinked

Author Summary
It is expected that genome-wide interaction analysis can be a possible source of finding heritability unexplained by current GWAS. However, the existing statistics for testing interaction have low power for genome-wide interaction analysis. To meet challenges raised by genome-wide interactional analysis, we develop a novel statistic for testing interaction between two loci (either linked or unlinked) and validate the null distribution and the type I error rates of the new statistic through simulations. By extensive power studies we show that the developed novel statistic has much higher power to detect interaction than the classical logistic regression. To provide evidence of gene-gene interactions as a possible source of the missing heritability unexplained by the current GWAS, we performed the genome-wide interaction analysis of psoriasis in two independent studies. The preliminary results identified 44 and 211 pairs of SNPs showing significant evidence of interactions with FDR,0.001 and 0.001,FDR,0.003, respectively, which were common in two independent studies. These included five interacting pairs of SNPs, some of which were located in the target sites: LST1/NCR3, CXCR5/BCL9L and GLS2 of miR-324-3p, miR-433, and miR-382, and 15 pairs of interacting SNPs that had nonsynonymous substitutions.
were considered. The statistics for testing interaction are usually motivated by the measure of interaction. The widely used logistic regression methods for detection of gene-gene interaction are based on then odds-ratio measure of interaction. Traditional additive and multiplicative odds ratio measures of interaction are defined in terms of genotypes at two loci. In this report, a novel statistic for testing interaction between two loci is based on multiplicative odds-ratio measures defined in terms of pseudohaplotypes. For the convenience of presentation, we first briefly introduce the odds ratio interaction measure in terms of genotypes, alleles, and then present the odds ratio measure in terms of pseudohaplotypes.

Genotype-Based Odds Ratio Multiplicative Interaction Measure
Consider two loci: G and H. Assume that the codes G~1(G~0) and H~1(H~0) denote whether an individual is a carrier (non-carrier) of the susceptible genotypes at the loci G and H, respectively. Let D denote disease status where D~1(D~0) indicates an affected (unaffected) individual. Consider the following logistic model: The odds-ratio associated with G for nonsusceptible genotype at the locus H (H~0) is defined as The odds-ratio associated with susceptibility at G and H compared to the baseline category G~0 and H~0 is then computed as OR GH~P (D~1DG~1,H~1)=P(D~0DG~1,H~1) P(D~1DG~0,H~0)=P(D~0DG~0,H~0) : The odds for baseline category G~0 and H~0 are determined as OR b~P (D~1DG~0,H~0) P(D~0DG~0,H~0) : From equation (1) If OR GH~O R G OR H , i.e., there is no interaction between loci G and H, then I GH~0 . This shows that the logistic regression coefficient for interaction term b GH is equivalent to the interaction measure defined as log odds-ratio. The interaction measure I GH can also be written as I GH~l og P(G~1,H~1jD~1)P(G~0,H~0jD~1) P(G~1,H~0jD~1)P(G~0,H~1jD~1) { log P(G~1,H~1jD~0)P(G~0,H~0jD~0) P(G~1,H~0jD~0)P(G~0,H~1jD~0) : The values of odds-ratio defined in terms of genotypes depends on how to code indicator variables G and H. Suppose that alleles G 1 and H 1 are alleles that increase disease risk. For a recessive model, G is coded as 1 if the genotype is G 1 G 1 , otherwise, G is coded as 0. For a dominant model, G is coded as 1 if the genotypes are either G 1 G 1 or G 1 G 2 , otherwise G is coded as 0. The indicator variable H can be similarly coded. However, in real data analysis, the disease models are unknown. Especially, the types of two-locus disease models are large [25]. We may have a large number of possible coding, and many of them may have larger numbers of degrees of freedom than the allelic model.

Allele-Based Odds Ratio Multiplicative Interaction Measure
Similar to the odds ratio for genotypes, we can define odds-ratio in terms of alleles. Let P(D~1DG i ,H j ) be the probability that an individual becomes affected given they have genotype G i =G k at locus G and H j =H I at locus H, where G k is either G 1 or G 2 (i.e. G k is a member of the set {G 1 ,G 2 }) and H I is either H 1 or H 2 (i.e. H I is a member of the set {H 1 ,H 2 }). We can similarly define P(D~0DG i ,H j ). We then can determine the odds-ratio associated with the allele G 1 at the G locus and allele H 1 at the H locus compared to the baseline G 2 =H 2 as Similarly, we measure the odds-ratio associated with the alleles G 1 =H 2 and G 2 =H 1 , respectively as : Similar to genotype, we can define a multiplicative interaction measure in terms of log odds-ratio for allele as which is equivalent to R~P (G 1 ,H 1 jD~1)P(G 2 ,H 2 jD~1) P(G 1 ,H 2 jD~1)P(G 2 ,H 1 jD~1) and S~P (G 1 ,H 1 jD~0)P(G 2 ,H 2 jD~0) P(G 1 ,H 2 jD~0)P(G 2 ,H 1 jD~0) : The ''fast-epistasis'' test statistic in PLINK (http://pngu.mgh. harvard.edu/,purcell/plink/index.shtml) is defined as where SE(R) and SE(S) denote the standard deviation of R and S, respectively. Absence of interaction is implied if and only if This is the basis of the ''fast-epistasis'' test in PLINK.

Haplotype-Based Odds Ratio Multiplicative Interaction Measure
Suppose that the locus G has two alleles G 1 and G 2 and the locus H has two alleles H 1 and H 2 . Let P A G1 ,P A G2 ,P A H1 ,P A H2 and P N G1 ,P N G2 ,P N H1 ,P N H2 be the frequencies of the alleles G 1 ,G 2 ,H 1 ,H 2 in the cases and controls, respectively. For the discussion of convenience, we introduce a terminology of ''pseudohaplotype''. When two loci are linked, a pseudohaplotype is defined as the regular haplotype. When two loci are unlinked, a pseudohaplotype is defined as a set of alleles that are located in the same paternal or maternal chromosomes. The frequencies of a pseudohaplotype can be estimated by the classical methods for estimation of haplotype frequencies such as Expectation Maximization (EM) Algorithms. For simplicity, hereafter we will not make distinction between the haplotype and pseudohaplotype. When two loci are unlinked, a haplotype is understood as a pseudohaplotype. Let P A 11 ,P A 12 ,P A 21 , P A 22 and P N 11 ,P N 12 ,P N 21 , P N 22 denote the frequencies of haplotypes G 1 H 1 ,G 1 H 2 ,G 2 H 1 and G 2 H 2 in the cases and controls, respectively. We define a penetrance of the haplotype G i H j as the probability that an individual becomes affected given they have phased genotype G i H j =G k H l ,k~1,2,l~1,2. Let f ijkl be the penetrance of an individual with the genotype G i H j =G k H l , h 11 ,h 12 ,h 21 and h 22 be the penetrance of the haplotypes G 1 H 1 ,G 1 H 2 ,G 2 H 1 and G 2 H 2 , respectively. The penetrance of the haplotype G i H j can be mathematically defined as where P 11 ,P 12 ,P 21 and P 22 are the population frequencies of the haplotypes G 1 H 1 ,G 1 H 2 ,G 2 H 1 and G 2 H 2 , respectively.
G~i and H~j represent a genotype coding scheme. Their represented genotypes depend on the specific genotype coding scheme. It should be noted that the haplotype G i H j and G~i and H~j have different meanings. By the same idea in defining genotype-based odds ratio in terms of penetrance of combinations of genotypes, we can determine the odds-ratio associated with the haplotypes G 1 H 1 compared to the baseline haplotype G 2 H 2 in terms of penetrance of the haplotypes as Similarly, we calculate the odds-ratio associated with the haplotypes G 1 H 2 and G 2 H 1 , respectively, as It is noted that replacing G~i and H~j in the definition of oddsratio in terms of genotypes by G i H j leads to the definition of oddsratio based on the haplotypes. However, the values and biological meanings of these two types of odds-ratios are different. Similar to genotypes, we can compute a multiplicative interaction measure in terms of log odds-ratio for haplotypes as In the absence of interaction, we have The multiplicative odds-ratio interaction measure in equation (3) is defined by the penetrance of the haplotypes. From case-control data it is difficult to calculate the penetrance of the haplotypes. However, we can show that the multiplicative odds-ratio interaction measure in equation (3) can be reduced to (Text S1, Appendix A) There are many algorithms and software to infer the haplotype frequencies in cases and controls. Therefore, we can easily calculate the multiplicative odds-ratio interaction measure by equation (4). It can be seen from equation (4)  in the controls are equal.
To gain understanding the multiplicative odds-ratio interaction measure, we study several special cases.
Case 1. One of two loci is a marker. If we assume that the locus H is a marker and is not associated with disease, then we have which implies that

21
: Thus, we obtain I H GH~0 . In other words, if the locus H is a marker, there is no interaction between two loci G and H. The interaction measure I H GH between two loci should be equal to zero. Hence, our multiplicative odds-ratio interaction measure correctly characterizes the marker case.
Case 2. Logistic regression interpretation. We define two indicator variables: Then four haplotypes at two loci can be coded as follows: It follows from the logistic regression model in equation (1) that where odds-ratios OR G and OR H are defined in terms of alleles, i.e.
OR G~P (D~1jG 1 )=P(D~0jG 1 ) P(D~1jG 2 )=P(D~0jG 2 ) and OR H~P (D~1jH 1 )=P(D~0jH 1 ) P(D~1jH 2 )=P(D~0jH 2 ) Therefore, the haplotype multiplicative odds-ratio interaction measure I H GH is equal to I H GH~bGH , which has the same form as that in equation (2-B). This indicates that if the coding for the genotypes in the genotype multiplicative odds-ratio interaction measure I GH is replaced by the coding for the haplotypes in equation (5) then we can obtain the haplotype multiplicative oddsratio interaction measure.

Test Statistics
In the previous section we defined the haplotype multiplicative odds-ratio interaction measure, which can be estimated by haplotype frequencies in cases and controls. By the delta method, we can obtain the variance of the estimator of the haplotype oddsratio interaction measure [26]: where n A and n G are the number of sampled individuals in cases and controls. By the standard asymptotic theory we can define the haplotype odds-ratio interaction measure-based statistic for testing interaction between two loci: whereP P A 11 ,P P A 12 ,P P A 21 ,P P A 22 andP P N 11 ,P P N 12 ,P P N 21 ,P P N 22 are the estimators of the corresponding haplotype frequencies in cases and controls, respectively. When sample sizes are large enough to ensure application of large sample theory, T IH is asymptotically distributed as a central x 2 (1) distribution under the null hypothesis of no interaction between two loci. Under an alternative hypothesis of of interaction between two loci being present, the statistic T IH is asymptotically distributed as a noncentral x 2 (1) distribution with noncentrality parameter proportional to the haplotype multiplicative odds-ratio interaction measure. This statistic can be applied to both linked and unlinked loci. As we explained in Text S1, Appendix B, the proposed statistic T IH is different from the ''fast-epistasis'' test in PLINK.
For the unlinked loci, we can use case only design [27,28] to study interaction between two loci in which equation is reduced to Results

Null Distribution of Test Statistics
In the previous sections, we have shown that when the sample size is large enough to apply large sample theory, the distribution of the statistic T IH for testing the interaction between two loci under the null hypothesis of no interaction between them is asymptotically a central x 2 (1) distribution. To examine the validity of this statement, we performed a series of simulation studies. MATLAB was used to generate two-locus genotype data of the sample individuals. A total of 100,000 individuals from a general population with an allele frequency P(G 1 )~0:4, P(H 1 )~0:3, haplotype frequency P(G 1 H 1 )~0:1 and disequilibrium coefficient d~P(G 1 H 1 ){P(G 1 )P(H 1 )~{0:02 were generated. A total of 10,000 simulations were repeated. Type I error rates were calculated by random sampling 500-1,000 individuals as cases and controls from the general population. Table 1 and Table 2 show that the estimated type I error rates of the statistic T IH for testing interaction between two loci, assuming OR G~O R H~1 and OR G~O R H~2 , were not appreciably different from the nominal levels a~0:05, a~0:01 and a~0:001. To further examine the validity of the test statistic, we constructed Quantile-quantile (Q-Q) plots of the test statistic in datasets 1 and 2 shown in Figures 1A and 1B, where the P-values of the tests were plotted (as 2log10 values) as a function of p values from the expected null distribution. Since the total number of all possible pair-wise tests for interaction between SNPs is too large to store all the results in computer we only stored P-values ,1: Consequently, Q-Q plots started with 4. Figures 1A and 1B showed good agreement with the null distribution.

Power Evaluation
To evaluate the performance of the statistic T IH for detection of interaction between two loci, we compared the power of the statistic T IH to that of the logistic model and the ''fast epistasis'' test in PLINK. Power was calculated by simulation. A total of 1,000,000 individuals from a general population with allele frequencies P(G 1 )~0:2, P(H 1 )~0:3 and P(G 1 H 1 )~0:1 and disequilibrium coefficient d~P(G 1 H 1 ){P(G 1 )P(H 1 )~0:04 were generated. Two-locus disease models were used to generate cases and controls, and summarized in Table 3 where odds-ratio was defined in terms of genotypes. We considered three types of genotype coding. For a recessive model, homozygous wild type, heterozygous, and homozygous risk increasing genotypes were coded as 0, 0, 1, respectively. For a dominant model, homozygous wild type, heterozygous, and homozygous risk increasing genotypes were coded as 0, 1, and 1, respectively. For an additive model, they were coded as 0, 1, and 2, respectively. The genotype coding for the logistic regression matched the simulation model. The statistic T IH in equation (6) for the case-control version was used to evaluate the power. In the power simulations, we also assumed that OR G~1 and OR H~1 . An individual who is randomly sampled from the general population was assigned to case or control status depending on the two-locus disease models in Table 3. The process was repeated until a sample of 1,000 cases and 1,000 controls for the dominant and additive models, or a sample of 2,000 cases and 2,000 controls for the recessive model was obtained. A total of 10,000 simulations were repeated. In Figures 2A-2C, power comparisons among the logistic regression model, the ''fast-epistasis'' in PLINK and the statistic T IH under two-locus recessive|recessive disease model for significance levels a~0:05, a~0:01 and a~0:001, respectively are presented. In Figures 3A-3C, power comparisons among the logistic regression model, the ''fast-epistasis'' in PLINK and the statistic T IH under two-locus dominant|dominant disease model for significance levels a~0:05, a~0:01 and a~0:001, respectively are shown. In Figures 4A-4C, power comparisons between the logistic regression model and the statistic T IH under two-locus additive|additive disease model for significance levels a~0:05, a~0:01 and a~0:001, respectively are demonstrated. Several remarkable features emerge from these Figures. First, these power Figures indeed demonstrate that the power increases as the measure of the interaction between two loci increases. The power curves were plotted as a function of the traditional genotype odds ratio OR GH . We observed that the power curves were a monotonic increasing function of the genotype odds ratio OR GH . Therefore, the test statistic T IH can detect the strength of the interaction between two loci. Second, the test statistic T IH had much higher power to detect interaction between two loci than the logistic regression and the ''fast-epistasis'' test in PLINK. Third, the more complex the disease models were, the larger the differences in power between the test statistic T IH , the ''fast-epistasis'' test in PLINK and logistic regression that were observed.
When two loci are unlinked where we do not observe the allelic association between two loci in the population as a whole, our results also hold. We assumed the following allele and haplotype frequencies in the population: P(G 1 )~0:2, P(H 1 )~0:3 and P(G 1 H 1 )~0:06. Other parameters were defined as before. A total of 10,000 simulations were repeated to simulate the power of three statistics under three disease models with the significance level a~0:001. Figures 5A, 5B and 5C showed the power of three statistics for testing interaction between two unlinked loci under two-locus recessive|recessive, dominant|dominant, and additive |additive disease models, respectively. These Figures again demonstrated that the power of the test statistic T IH was still much higher than that of the logistic regression and the ''fast-epistasis'' test in PLINK. The conclusions still hold for the significance levels a~0:05 and a~0:01 (Data were not shown).

Application to Pathway-Based Genome-Wide Interaction Analysis of Psoriasis
To evaluate its performance for detection of interaction between two loci, the proposed test statistic T IH was applied to interaction analysis of two independent GWAS datasets of psoriasis which were downloaded from dbGaP. Psoriasis is a common chronic inflammatory skin disease affecting 2%-3% of the world population. Originally, the first study included 955 individuals with psoriasis and 693 controls, which is considered as dataset 1. The second replication study included 466 individuals with psoriasis and 732 controls, which is designated dataset 2. All cases and controls are of European origin [29][30][31]. After using PLINK [20] to check for contamination, cryptic family relationship and non-Caucasian ancestry, 123 samples were excluded. Subsequently we retained for analysis 915 cases and 675 controls from the first study and 431 cases and 702 controls from the second study. All 2,723 samples had been genotyped with the Perlegen 500K array. In the initial dataset, 451,724 SNPs passed quality control (call Table 2. Type I error rates of the statistic T IH to test for interaction between two loci, assuming OR G~O R H~2 . Since testing for all possible two-locus interactions across the genome in genome-wide interaction analysis requires extremely large computation, we conducted pathway-based genome-wide interaction analysis. We assembled 501 pathways from KEGG [24] and Biocarta (http://www.biocarta.com). The assignment of SNPs to a gene was obtained from NCBI human9606 database (version b129). We used the statistic T IH to test interactions of all possible pairs of SNPs located in genes within the assembled 501 pathways. The total number of SNPs in dataset 1 and dataset 2 being tested was 50,327 and 49,876, respectively. The serious problem in genome-wide interaction analysis is multiple testing. We used two strategies to tackle this problem. One is to use false discovery rate (FDR) [32] to declare significance of interaction. Another is replication of the findings in two independent studies, which enhances confidence in interaction tests [22]. We looked for consistent results across the two independent studies.
In total, 44 pairs of SNPs showed significant evidence of interactions with FDR,0.001, which roughly corresponds to the P-value ,1:0|10 {7 , in two independent studies (Table S1). These 44 pairs of SNPs were derived from 71 distinct SNPs located in 60 genes, including HLA-C, HLA-DRA, HLA-DPA1, LST1, MICB and NOTCH4. Of 44 pairs of SNPs, only one pair of interacting SNPs: rs2395471 and rs2853950 showed significant marginal association in two independent studies. An additional 211 pairs of SNPs with FDR less than 0.003 in the two studies is listed in Table  S2. These interacting SNPs were mainly located in 19 pathways, including a number of signaling pathways, and immune-related antigen processing and presentation as well as natural killer cell mediated cytotoxicity pathways ( Figure 6). Several remarkable features emerge from these results. First, although we can observe a few interactions between SNPs within a gene, the majority of  interactions occurred between genes that are often in different pathways. Since the number of SNPs typed within each gene was limited, it is unknown whether this is a general rule or just a special case. Second, a SNP in one gene might interact with multiple SNPs in multiple genes. For example, SNP rs3131636 in the gene MICB interacting with the SNPs rs915895, rs443198, rs3134929 in the gene NOTCH4, the SNP rs1052248 in the gene LAST1/ Natural cytotoxicity triggering receptor 3 (NCR3) and the SNP rs1799964 in the gene LTA/TNF. SNP rs1799964 in the gene LTA/TNF interacting with SNPs rs3131636, rs3132468 in the gene MICB, SNPs rs9268658 and rs3135392 in the gene HLA-DRA, SNP rs2227956 in the gene HSPA1L. However, this does not imply that multiple causal SNPs within a gene will interact with multiple causal SNPs within another gene. It is quite likely that this is due to LD between the SNPs within a gene. Third, although interacting SNPs did not form large connected networks, the interacting SNPs connected pathways into a large complicated network. This may imply that many genes and pathways are involved in the development of psoriasis. Fourth, upstream of many pathways included genes with interacting SNPs. For example,  genes MICB, CHRM3, HLA-DRA and CIITA, EPHB1 and EPHB2, LAMA1 and LANA5, ITGA1, LTBP1, TNF, and FGF20 that contain interacting SNPs are in the upstream of natural killer cell mediated cytotoxicity, calcium signaling pathway, antigen processing and presentation, axon guidance, ECM-receptor interaction pathway, focal adhesion, TGFB pathway, MAPK pathway and regulation of acting cytoskeleton, respectively. Fifth, most interacting SNPs are in introns and accounted for 77% of total interacting SNPs. Table 4 listed 15 pairs of interacting SNPs that have nonsynonymous substitutions. It is unknown how these nonsynon-ymous mutations are involved in the pathogenesis of psoriasis. From the literature we know that Plexin C1 receptor is a tumor suppressor gene for melanoma [33], NOTCH4 is involved in schizophrenia [34], Phosphodiesterase 4D (PDE4D) is associated with ischemic stroke [35], HLA-DRA is one of the HLA class II alpha chain genes that plays a central role in antigen processing, and neuregulin 1 (NRG1) has been implicated in diseases such as cancer, schizophrenia and bipolar disorder [36]. Table 5 includes five interacting pairs of SNPs, one of which falls in the microRNA (MiRNA) binding region. miRNAs, which are 22 nucleotide small RNAs and regulate gene expressions by Figure 5. Power of the statistics for testing interaction between two unlinked loci. (A) The power of the test statistic T IH , the ''fastepistasis'' in PLINK and logistic regression analysis for testing interaction between two unlinked loci as a function of traditional odds-ratio OR GH under a two-locus recessive|recessive disease model, where the number of individuals in both the case and control groups is 2,000, the significance level is 0.001, and the odds-ratios at two loci were OR G~O R H~1 . (B) The power of the test statistic T IH , the ''fast-epistasis'' in PLINK and logistic regression analysis for testing interaction between two unlinked loci as a function of traditional odds-ratio OR GH under a two-locus dominant|dominant disease model, where the number of individuals in both the case and control groups is 1,000, the significance level is 0.001, and the odds-ratios at two loci were OR G~O R H~1 . (C) The power of the test statistic T IH , the ''fast-epistasis'' in PLINK and logistic regression analysis for testing interaction between two unlinked loci as a function of traditional odds-ratio OR GH under a two-locus additive|additive disease model, where the number of individuals in both the case and control groups is 1,000, the significance level is 0.001, and the odds-ratios at two loci were OR G~O R H~1 . doi:10.1371/journal.pgen.1001131.g005 The power of the test statistic T IH , the ''fast-epistasis'' in PLINK and logistic regression for testing interaction between two linked loci analysis as a function of traditional odds-ratio OR GH under a two-locus additive|additive disease model, where the number of individuals in both the case and control groups is 1,000, the significance level is 0.05, and the odds-ratios at two loci were OR G~O R H~1 . (B) The power of the test statistic T IH , the ''fast-epistasis'' in PLINK and logistic regression for testing interaction between two linked loci analysis as a function of traditional odds-ratio OR GH under a two-locus additive|additive disease model, where the number of individuals in both the case and control groups is 1,000, the significance level is 0.01, and the odds-ratios at two loci were OR G~O R H~1 . (C) The power of the test statistic T IH , the ''fast-epistasis'' in PLINK and logistic regression for testing interaction between two linked loci analysis as a function of traditional odds-ratio OR GH under a two-locus additive|additive disease model, where the number of individuals in both the case and control groups is 1,000, the significance level is 0.001, and the odds-ratios at two loci were OR G~O R H~1 . doi:10.1371/journal.pgen.1001131.g004 pairing the miRNA seed region with the target sites, have been implicated in many biological processes including the immune response, biogenesis and tumorigenesis [37]. Mutations in the target sites will affect miRNA activity. A number of studies have identified polymorphisms in the miRNA target sites associated with the diseases [37]. Interestingly, we identified four SNPs in the miRNA (miR-324-3p, miR-433, and miR-382) target sites which interact with five SNPs to contribute to psoriasis. In previous studies, miR-382 has been associated with dermatomyositis, Duchenne muscular dystrophy and Miyoshi myopathy [38], miR-433 and miR-324 with lupus nephritis [39] and miR-433 with Parkinson's disease [40].
Some researchers suggest that in genome-wide interaction analysis only SNPs with large or mild marginal genetic effects should be tested for interaction. To examine whether this strategy will miss detection of interacting SNPs, we showed in Table 6 the 20 top pairs of interacting SNPs and in Table S3 all pairs of interacting SNPs with FDR less than 0.003. Surprisingly, 75% of SNPs with P-values (in dataset 1) larger than 0.2 and 44% of SNPs with P-values larger than 0.5 in two studies were observed in Table  S3. Table 6 and Table S3 strongly demonstrated that while both SNPs did not demonstrate significant evidence of marginal association, they did show significant evidence of interaction.
To further evaluate the performance of the proposed statistic T IH , in Table 7 and Table S4 we list P-values for testing interaction calculated by the statistic T IH , the ''fast-epistasis'' in PLINK and logistic regression using genotype coding. In Table 7 the 20 top pairs of interacting SNPs and in Table S4 the results of 233 pairs of interacting SNPs are presented. The P-values for interaction calculated by the statistic T IH are much smaller than those from the ''fast-epistasis'' in PLINK and the logistic regression using genotype coding (Table 7 and Table S4). Moreover, the ''fast-epistasis'' in PLINK and the logistic regression coded by genotype detect very few interactions that can be replicated in two independent studies (Table 7 and Table S4). In fact, our results for all tested SNPs in 501 pathways showed that the ''fast-epistasis'' in PLINK and logistic regression coded by genotypes detected very few interactions that can be replicated in two studies (data not shown).
Eighteen significantly interacting SNPs identified by Bonferroni correction were listed in Table 8. In dataset1, the total number of SNPs for testing interaction was 50,327. The P-values for declaring interaction between SNPs after Bonferroni correction was 3:95|10 {11 . We found that there were 2,210 significant interactions with P-values less than 3:51|10 {11 in the dataset 1. Then, interaction for all these 2,210 pairs of SNPs in the dataset 2 Figure 6. Interacting SNPs that were located in 19 pathways formed a network. Each pathway was represented by an ellipse with the number. The SNPs were represented by nodes and placed insight their located pathways. Nearby each SNP there was its RS number and the name of its located gene. The pathway and its harbored SNPs were labeled by the same color. The interacting SNPs were connected by the solid light green lines. doi:10.1371/journal.pgen.1001131.g006 was examined. The P-values for declaring interaction between SNPs after Bonferroni correction in dataset 2 was 2:26|10 {5 . We identified eight significant interactions that were replicated in dataset 2. Similarly, if we started with dataset 2, the total number of SNPs for testing interaction was 49,876. The P-values for declaring interaction between SNPs after Bonferroni correction was 4:02|10 {11 . Significant interactions with the P-values less than 3:87|10 {11 in dataset 2 were seen between 1,913 pairs of SNPs. Then, we tested for interaction for all these 1,913 pairs of SNPs in the dataset 1. The P-values for declaring interaction between SNPs after Bonferroni correction in dataset 1 was 2:61|10 {5 , and 10 significant interactions were detected that were replicated in the dataset 1. A total of 9 interactions were common in Table 8 and Table S1 and Table S2.

Discussion
The development of most diseases is a dynamic process of genegene and gene-environment interactions within a complex biological system. We expect that genome-wide interaction analysis will provide a possible source of finding missing heritability unexplained by current GWAS that test association individually. But, in practice, very few genome-wide interaction analyses have been conducted and few significant interaction results have been reported. Our aim is to develop statistical methods and computational algorithms for genome-wide interaction analysis which can be implemented in practice and provide evidence of gene-gene interaction. The purpose of this report is to address several issues to achieve this goal.
The first issue is how to define and measure interaction. Oddsratio is a widely used measure of interaction for case-control design. The odds-ratio based measure of interaction between two loci is often defined as a departure from additive or multiplicative odds-ratios of both loci defined by genotypes. The genotype-based odds-ratio does not explore allelic association information between two loci generated by interaction between them in the cases. Any statistics that are based on genotype defined odds-ratio will often have low power to detect interaction. To overcome this limitation, we extended genotype definition of odds-ratio to haplotypes and revealed relationships between haplotype-defined odds-ratio and haplotype formulation of logistic regression. To further examine the validity of this concept, we studied the distribution of the test statistic under the null hypothesis of no interaction between two either linked or unlinked loci. Through extensive simulation  (assuming allelic association in the controls), we show that the distribution of the haplotype odds-ratio-based statistic is close to a central x 2 (1) distribution even for small sample size and that type I error rates were close to the nominal significance levels.
The second issue is the power of the test statistic for genomewide interaction analysis. The genome-wide interaction analysis requires testing billions of pairs of SNPs for interactions. The Pvalues for ensuring genome-wide significance level should be very small. Therefore, developing statistics with high power to detect interaction is an essential issue for the success of genome-wide interaction analysis. As an alternative to the logistic regression and the ''fast-epistasis'' in PLINK, we presented a haplotype oddsratio-based statistic for detection of interaction between two loci and illustrated its power by extensive simulations. The power of the haplotype odds-ratio-based statistic ended up being a function of the measure of interaction and had much higher power to detect interaction than the ''fast-epistasis'' in PLINK and logistic regression.
The third issue is whether the interactions exist with no marginal association and how often they might occur in practice. Our data demonstrated that the majority of the significantly interacting SNPs showed no marginal association. Surprisingly, 75% of interacting SNPs with P-values (for testing marginal association) larger than 0.2 and 44% of interacting SNPs with Pvalues (for testing marginal association) larger than 0.5 in two studies were observed in our analysis. This strongly suggested that testing interaction for only SNPs with strong or mild marginal association will miss the majority of interactions.
The fourth issue is that of replication of the results. Genomewide interaction analysis involves testing billions of pairs of SNPs.
Even if after correction of multiple tests, the false positive results might be still high. To increase confidence in interaction test results, replication of interaction findings in independent studies is often sought. To date, very few results of genome-wide interaction analysis have been replicated. This begs the question whether the significant interaction can be replicated in independent studies. In this report, we show that interaction findings can be replicated in two independent studies.
The fifth issue is correction for multiple testing. Genome-wide interaction analysis often involves billions of tests, which would require an extremely small Bonferroni-corrected P-value to ensure a genome-wide significance level of 0.05. Replication of finding at such small P-values in independent studies is often extremely difficult. However, Bonferroni correction assumes that the tests are independent, yet many interaction tests are highly correlated. Correlations in the interaction tests come from two levels [23]. First, two pairs of SNPs may share a common SNP. Second, SNPs in the interaction tests may be dependent due to allelic association. The Bonferroni correction assuming independent tests will be overly conservative due to high correlations among the interaction tests. In this report, two strategies were used to tackle the multiple testing issues. The first is to use FDR to control type I error. The second is to replicate interaction finding. Replication allows us to detect the interactions that are frequent and consistent [22]. This approach still has the limitation that we still make independent assumption of the tests in calculation of FDR. Recently, Emily et al. (2009) [23] proposed to develop a Bonferroni-like correction for multiple tests based on the concept of the effective number of SNP pairs. The concept of the effective number of tests takes correlation among the tests into account and can be applied to both P-value  and FDR correction [41]. This may be a promising approach to multiple test corrections in the genome-wide interaction analysis. Although our data show that interactions can partially find the heritability of complex diseases missed by the current GWAS, they are still preliminary. Due to extremely intensive computations demanded by genome-wide interaction analysis we only tested interactions of a small set of SNPs which were located in the genes of 501 assembled pathways in a PC computer. The truly whole genome interaction analysis in which we will test for interactions between all possible pairs of SNPs across the genome has not been conducted. Gene-gene interaction is an important, though complex concept. The statistical interactions are scale dependent. There are a number of ways to define gene-gene interaction. How to define gene-gene interaction and develop efficient statistical methods and computational algorithms for genome-wide interaction analysis are still great challenges facing us. The main purpose of this report is to stimulate discussion about what are the optimal strategies for genome-wide interaction analysis. We expect that in coming years, genome-wide interaction analysis will be one of major tasks in searching for remaining heritability unexplained by the current GWAS approach.