A hidden two-locus disease association pattern in genome-wide association studies

Background Recent association analyses in genome-wide association studies (GWAS) mainly focus on single-locus association tests (marginal tests) and two-locus interaction detections. These analysis methods have provided strong evidence of associations between genetics variances and complex diseases. However, there exists a type of association pattern, which often occurs within local regions in the genome and is unlikely to be detected by either marginal tests or interaction tests. This association pattern involves a group of correlated single-nucleotide polymorphisms (SNPs). The correlation among SNPs can lead to weak marginal effects and the interaction does not play a role in this association pattern. This phenomenon is due to the existence of unfaithfulness: the marginal effects of correlated SNPs do not express their significant joint effects faithfully due to the correlation cancelation. Results In this paper, we develop a computational method to detect this association pattern masked by unfaithfulness. We have applied our method to analyze seven data sets from the Wellcome Trust Case Control Consortium (WTCCC). The analysis for each data set takes about one week to finish the examination of all pairs of SNPs. Based on the empirical result of these real data, we show that this type of association masked by unfaithfulness widely exists in GWAS. Conclusions These newly identified associations enrich the discoveries of GWAS, which may provide new insights both in the analysis of tagSNPs and in the experiment design of GWAS. Since these associations may be easily missed by existing analysis tools, we can only connect some of them to publicly available findings from other association studies. As independent data set is limited at this moment, we also have difficulties to replicate these findings. More biological implications need further investigation. Availability The software is freely available at http://bioinformatics.ust.hk/hidden_pattern_finder.zip.


Background
The development of DNA microchip technology has allowed the analysis of single nucleotide polymorphism (SNPs) on a genome-wide scale to identify genetic variants associated with diseases. Researchers have proposed many methods to investigate association patterns of complex diseases. Two recent reviews [1,2] presented detailed analyses on many popular methods and tools, such as multifactor dimensionality reduction (MDR) [3], Random Jungle [4], Bayesian epistasis association mapping (BEAM) [5] and PLINK [6]. MDR is a popular non-parametric approach for detecting all possible k-way combinations of SNPs that interact to influence disease traits. Random Jungle (i.e., Random Forest [7]), is to solve classification and regression problems. In random forest, decision trees are combined to produce accurate predication. Its ability to handle the high dimensional problems in GWAS has been shown in [8,9]. BEAM designs a Bayesian marker partition model which classifies SNP markers into three types: SNPs unassociated with the disease, SNPs contributing to the disease susceptibility independently, and SNPs influencing the disease risk jointly with each other. In this model, a first order Markov chain is designed for the accommodation of correlation between adjacent SNPs. Markov Chain Monte Carlo (MCMC) sampling is used to optimize the posterior probability of the model. In addition, the "B-statistic" designed in BEAM can be used in the frequentist hypothesis-testing framework. PLINK provides a toolkit for flexible analyses, in which various statistical tests for single-locus analysis, haplotype analysis and allelic-based interaction analysis are implemented. Recently, a new method named "BOOST" [10] allows examination of all pairwise interactions in genome-wide case-control studies. As a result, many genetic susceptibility determinants have been mapped.
However, there is another type of association pattern, which often occurs within local regions in the genome and may not be detected by these methods. This association pattern involves multiple correlated SNPs with neither strong marginal effects nor strong interaction effects. But they can jointly display strong associations. Here we use some simple regression models to explain this association pattern. Suppose we have two dependent variables, X 1 and X 2 , and one independent variable Y. We can fit two regression models (or logistic regressions for case-control data), Y ∼β 1 X 1 and Y ∼β 2 X 2 , to test the association significance of these two variables. Hereβ 1 andβ 2 are named as marginal coefficients. If these two marginal coefficients are very small, single variable analysis methods will consider them statistically insignificant and ignore them.
However, if X 1 correlates with X 2 , fitting the model Ỹ b 1 X 1 + b 2 X 2 may identify a new association pattern with b 1 and b 2 (named as bivariate regression coefficients) being significantly larger thanβ 1 andβ 2 . This phenomenon is referred to as unfaithfulness. It means that the marginal effects of correlated variables do not express their significant joint effects faithfully due to the correlation cancelation [11]. Figure 1 provides some synthetic examples to show the unfaithfulness involving two variables. There are four scenarios to illustrate the relationship between marginal coefficients (marked using red color) and bivariate regression coefficients. The first scenario ( Figure 1:(a)) is a reference case that involves no correlations between X 1 and X 2 . The marginal coefficientsβ 1 andβ 2 are equal to the bivariate regression coefficients b 1 and b 2 , respectively. In the second scenario ( Figure 1:(b)), X 1 is positively correlated with X 2 . The marginal coefficients are bigger than the bivariate regression coefficients. In the third scenario (Figure 1: (c)), X 1 is negatively correlated with X 2 . The marginal coefficientβ 1 andβ 2 could be significantly smaller than the bivariate regression coefficients b 1 and. b 2 In the the fourth scenario ( Figure 1:(d)), X 1 is positively correlated with X 2 . But the sign of b 1 is the opposite of the sign of β 2 . The correlation effect in the third scenario and the fourth scenario causes the unfaithfulness. In mathematics, the relationship between the marginal coefficients and the bivariate regression coefficients is formulated as where E( β 1 ) is the expectation of the marginal coefficientβ 1 , r(X 1 , X 2 ) is the population correlation between Figure 1 Illustration of unfaithfulness in association studies. There are three regression models in each scenario: Y~b 1 X 1 + b 2 X 2 , Y ∼β 1 X 1 and Y ∼β 2 X 2 . In this figure, the marginal coefficientβ 1 andβ 2 are shown as projections (marked with bold red color) of Y on X 1 and X 2 , respectively. (a) X 1 is not correlated with X 2 . (b) X 1 is positively correlated with X 2 . (c) X 1 is negatively correlated with X 2 . (d) X 1 is positively correlated with X 2 but the sign of b 1 is the opposite of the sign of b 2 . Scenario (c) and Scenario (d) illustrate unfaithfulness. X 1 and X 2 . The marginal coefficients depend on their bivariate regression coefficients as well as the variable correlation, as we illustrated in Figure 1. We will give more explanation on this relationship in the high dimensional setting in the discussion section.
From the statistical point of view, different correlation patterns could cause the marginal coefficientsβ 1 andβ 2 significantly different from the bivariate regression coeffcients b 1 and b 2 (shown in Figure 1). In fact, the issue of unfaithfulness has been discussed in the causality literature [12]. In GWAS, the correlation among SNPs arises due to the linkage disequilibrium pattern of the genome. A natural question arises: Does the issue of unfaithfulness occur in GWAS?
To answer this question, a computational method for detecting associations masked by unfaithfulness is necessary. In this paper, we propose a simple method to detect such associations involving two SNPs. It can evaluate each SNP pair in genome-wide case-control studies in a fast manner. We have applied our method to analyze seven data sets from the Wellcome Trust Case Control Consortium (WTCCC). The experimental results show that these associations widely exist in GWAS. In this work, we only handle the unfaithfulness issue involving two SNPs while the unfaithfulness can exist among a large number of markers. The detection of associations involving three or more SNPs is too time-consuming and beyond the scope of this work.

Experiment on simulation study
The simulation study is designed to compare our proposed method with other three methods for detecting associations in the presence of unfaithfulness. These three methods include the marginal association test (single-locus analysis), Lasso [11,13] and BEAM [5]. The reasons that we choose these methods for comparison are as follows: • Marginal association test is used in almost every GWAS due to its simplicity and effectiveness.
• Lasso is a shrinkage and selection method for (generalized) linear regression. It imposes a sparsity constraint (i.e., only a small fraction of variables are relevant) and uses L 1 penalty to eliminate irrelevant variables. Fast algorithms are available for Lasso. Thus, it can simultaneously analyze a huge number of variables. It is very popular in genetics [11,[14][15][16].
• BEAM has the capability of detecting both marginal associations and interactions in large-scale data sets. It uses first order Markov chain to accommodate the correlation between adjacent SNPs.
The details about the parameter settings in simulation are provided in the method section. In our simulation study, we only handle the unfaithfulness involving two associated variables X 1 and X 2 by using β 1 > 0, b 2 < 0 and r(X 1 , X 2 ) > 0 as illustrated in Figure 1(d). The marginal coefficientsβ 1 andβ 2 will be small due to the cancelation given by Equation (1). When β 1 > 0, b 2 > 0 and r(X 1 , X 2 ) < 0 the unfaithfulness also happens. This corresponds to a situation that the minor alleles of both X 1 and X 2 increase the diseases risk but X 1 and X 2 are negatively correlated, as illustrated in Figure 1(c).
The results in Figure 2 indicate that it is difficult for existing methods to detect the association masked by unfaithfulness while our proposed method achieves reasonable performance. Specifically, the poor performance of the marginal association test is not surprising since the marginal effects are weak in the presence of unfaithfulness. Although Lasso can simultaneously analyze all SNPs, it still suffers from the difficulty of detecting associations masked by unfaithfulness. This agrees with the analysis result in [11]. BEAM has a better performance, which should be attributed to its first order Markov chain designed for the accommodation of correlation. But its performance is still not comparable with the performance of our proposed method in most settings.
Another interesting point is that the statistical power of existing methods decreases as the linkage disequilibrium (LD) r 2 increases. Although our proposed method also degrades its performance when LD increases, it maintains a relatively high power for strong LD (r 2 = 0.7).

Experiment on seven data sets from WTCCC
We have applied our method to analyze the data sets (14,000 cases in total and 3,000 shared controls) from the WTCCC [17]. WTCCC studies seven common human diseases, including bipolar disorder (BD), coronary artery disease (CAD), Crohn's disease (CD), hypertension (HT), rheumatoid arthritis (RA), type 1 diabetes (T1D) and type 2 diabetes (T2D). These data sets are generated using the affymetrix 500 K chip. We first apply a similar quality control procedure as suggested in [17] to pre-process the data. The numbers of remaining SNPs for seven data sets are around 360,000. In current stage, BEAM cannot directly handle these data sets [2]. Table 1 lists the numbers of identified two-locus associations masked by unfaithfulness under three statistical significance thresholds with and without the distance threshold for seven data sets. It shows that the unfaithfulness widely exists in these data sets. Some associations masked by unfaithfulness involve SNPs with at least 1 M base pair distance. However, all of them are located in the major histocompatibility complex (MHC) region (The MHC region encodes a large number of genes. It has extensive polymorphism and linkage disequilibrium with the long distance [18]). Therefore, the results in Table 1 provide the evidence that this association pattern typically occurs in local area. These results also suggest that using the local search can speed up the whole process in the future.
From the identified associations, we further conduct the gene mapping and identify some suspicious genes closely related with the disease traits. Table 2 and Table 3 report the unadjusted single-locus P-values, the unadjusted joint P-values, the marginal coefficients and joint bivariate coefficients for these associations. The other details are listed in the supplementary document (Additional file 1). These identified associations coincide with Figure 1(d). To date, we can only connect some identified associations to publicly available results from other association studies. Many identified association patterns still remain unexplained. In the following, we explain the details of some associations that are confirmed by other studies.

Bipolar disorder (BD)
Among associated SNP pairs identified from the BD data set, we find two suspicious SNP pairs, (rs668860, rs10873672) and (rs668860, rs6691970). The unadjusted P-values for these two SNP pairs are 4.885 × 10 -15 and 6.217 × 10 -15 , respectively. They are still significant after Bonferroni correction. However, none of these three SNPs (rs668860 is involved in both pairs) was reported in [17] because their marginal effects are too weak to be detected by the single-locus association test. The unadjusted P-values for these three SNPs based on the single-locus association test are 0.053, 0.245 and 0.216, respectively. All three SNPs reside in the intron of gene MCOLN2. The protein Mucolipin-2, encoded by gene MCOLN2 and also known as TRPML2 (transient receptor potential cation channel, mucolipin subfamily, member 2), has been confirmed to have strong associations with bipolar disorder in a family-based association study [19]. To our knowledge, this is the first identified association between the MCOLN2 gene and the bipolar disorder risk in a population-based association study. Figure 3 shows the joint distributions of the pair (rs668860, rs10873672) (The other pair shares a similar pattern.) in cases and controls and the corresponding odds ratios. The genotype combination "CT/TT" has a significantly higher odds ratio than other genotype combinations. Further investigations of the MCOLN2 gene may help identify the causes of bipolar disorder disease.
We further use BEAGLE [20] to impute the SNPs in this local area so that we can see the enriched signals after imputation. This region includes 300 SNPs. It begins with the SNP rs1030933 and ends with the SNP rs1837329. After imputation, we analyze the imputed data and the result is given in Figure 4. Figure 4(a) shows the enriched signal. The intensity shows -log 10 P values given by the joint regression (P-value is calculated based on χ 2 df =4 ). Figure 4(b) shows the LD structure of this local area. Figure 4(c) shows the -log 10 P values obtained using single-SNP analysis (P-values are calculated based χ 2 df =2 ). Figure 4(d) shows the locations of rs668860, rs10873672 and rs6691970. From Figure 4, we can see the following: • Although this region is in strong LD (see Figure 4 (b)), association masked by unfaithfulness does not happen across the entire area. This shows that this type of asssociation not only depends on the correlation structure but also depends on the effects of the SNPs, as we illustrated in Figure 1 (also see Equation  1).
• From Figure 4(c), the marginal effects of the imputed SNPs are very weak. This indicate that this type of association is not caused by some ungenotyped causative SNPs. Instead, it is a genuine effect.

Coronary artery disease (CAD)
We identify four suspicious associations involving five SNPs. The unadjusted P-values for these four associations range from 2.310 × 10 -13 to 5.551 × 10 -15 . The unadjusted single-locus P-values for five SNPs involved in these five associations indicate that they do not have noticeable marginal effects. All five SNPs reside in the intron of gene FSIP1 (fibrous sheath interacting protein 1). We have not found evidence to directly connect gene FSIP1 with the coronary artery disease. However, the LD analysis identifies a well studied gene THBS1 Dist -the physical distance threshold between two SNPs is at least 1 Mb. This threshold is used to see how many unfaithfulness associations involving two remote loci (thrombospondin-1), which is centromeric to gene FSIP1 and has been confirmed to increase the risk of coronary artery disease in many studies [21][22][23]. It would be of great interest to investigate gene FSIP1 in determining genetic susceptibility to coronary artery disease.
Here we also show the enriched signals obtained from the imputation. Figure 5(a) shows the -log 10 P given by the joint regression. Figure 5(b) shows the LD structure (r 2 ) in that region. Figure 5(c) shows the -log 10 P of single SNP analysis. Figure 5(d) shows the locations of the genotyped SNPs which are listed in Table 2. Again, the marginal effects of the imputed SNPs are weak. We see clearly that the signal of unfaithfulness appears in the block-like manner.

Type 1 diabetes (T1D)
Most identified associations from the T1D data set are linked with the MHC region. The MHC region at chromosomal position 6p21 encodes many genes (such as HLA-DQB1 and HLA-DRB1) that have been associated with type 1 diabetes [17,24] by using the singlelocus test. However, it is still unclear which and how many loci within the MHC region determine T1D susceptibility because of the functional complexity of this small human genome segment. The MHC region has been connected with more than 100 diseases, such as diabetes, rheumatoid arthritis, psoriasis, asthma and various autoimmune disorders. Our results provide additional information to locate disease-associated loci. Concretely, one suspicious association involves SNP rs1058318 and SNP rs2252745. The unadjusted P-value of this association is 1.326 × 10 -12 . The unadjusted single-locus P-values of rs1058318 and rs2252745 are 0.074 and 0.840, respectively. SNP rs1058318 resides in the intron region of gene GNL1 and SNP rs2252745 resides in the intron region of gene PPP1R10. Both genes are Table 3 Regression coefficients of those associations listed in Table 2 Disease SNP X pβ 1 (z Value) SNP X qβ 2 (z Value)  The distribution of genotypes of rs668860 and rs10873672 in control samples. Right panel: The estimated odds ratio for the combination of rs668860 and rs10873672. The odds ratio of genotype combination "AA/TT" is used as reference. The genotype combination "TT/CT" has significantly higher odds ratio than other genotype combinations. located in the MHC region and adjacent to each other. Gene GNL1 belongs to the HLA-E class. The locus in HLA-E has been strongly associated with type 1 diabetes [25]. The detailed examination of the relationship between gene GNL1 and gene PPP1R10 may provide some new insights in studying the causes of type 1 diabetes.

Hypertension (HT)
Among associations identified from the HT data set, we find one suspicious pair involving SNP rs2300390 and SNP rs12482676. The unadjusted P-value is 2.442 × 10 -15 . The unadjusted single-locus P-values for rs2300390 and rs12482676 are 0.460 and 0.061, respectively. Both SNPs reside in the intron of gene RCAN1. Gene RCAN1 mainly functions as a regulator of calcineurin. Calcineurin participates in many cellular and tissue functions. Its abnormal expression is associated with many diseases including hypertension [26].
Crohn's disease (CD), rheumatoid arthritis (RA) and type 2 diabetes (T2D) Currently, we have difficulties to connect the identified associations of CD, RA and T2D to publicly available findings from other association studies. Their biological implications need to be further explored.

Experiment on the Illumina data sets from other independent studies
We further analyze the Crohn's Disease data set [27], in which 308,332 autosomal SNPs were assayed on the Illumina HumanHap300 chip. After a standard quality control (the proportion of miss values ≤ 10%, the minor allele frequency ≥ 5% and the P-value of Hardy-Weinberg equilibrium ≥ 0.0001), the number of remaining SNPs is 291,964.
We apply our method to this data set and do not find any significant associations masked by unfaithfulness. Our explanation is that Illumina chip uses the tagSNP design and the correlation among SNPs is less than that of Affymetrix 500 K chip used by WTCCC. This result indicates that it is unlikely to detect associations masked by unfaithfulness using the tagSNP design.
In order to check if imputation helps in identifying significant association masked by unfaithfulness, we focus on the SNP regions in which we have identified associations from the WTCCC CD data set (Additional file 1: Table S3), impute the corresponding SNP data from [27], and re-run our analysis. Unfortunately, we fail to replicate those findings in the WTCCC CD data set. We have examined the imputation result carefully. At those local regions we are interested in, few SNPs are directly genotyped. In the hapmap data, hundreds of SNPs locate in those areas. This implies hundreds of SNPs need to been imputed using the information coming from the reference panel. In fact, the frequencies of those imputed haplotypes are almost the same in cases and controls. This is probably the reason that we cannot replicate those findings. Hopefully, next-generation sequencing will provide high resolution SNP data to resolve this issue. Another important reason may be that these two CD data sets are from different populations (one comes from Europe, another comes from north America). Similarly, we have analyzed another RA data set [28] from Genetic Analysis Workshop 16. This data set is acquired from North American population. The SNPs are genotyped by the Illumina chip. We also have difficulties to replicate the findings of the RA data set from WTCCC. We hope we can get access to more data sets to verify our results in the future.

Discussion
The unfaithfulness issue in the high dimensional feature space In the high dimensional feature space, many features could correlate with each other by chance, which leads to the existence of unfaithfulness and poses a great challenge on feature selection and association analysis. In this work, we only handle the unfaithfulness issue involving two variables (SNPs), while the unfaithfulness can exist among a huge number of variables. The relationship between the marginal coefficient (β i in Y ∼β i X i ) and the regression coefficient (β i in b 1 X 1 +···+b p X p +···+b s X s ) is given as follows [11]: where E( β p ) is the expectation of marginal coefficient, r(X q , X p ) is the population correlation between X q and X p . If β p ≈ − 1≤q≤s,q =p β q ρ(X q , X p ) , thenβ p can be close to zero no matter how large b p is. In addition, the number of involved variables could be very big. To exclude the effect of unfaithfulness in feature selection, Fan and Lv [29] had to make an assumption that there is a C > 0 such that |β p | ≥ C|β p | for p = 1,..., s, and then proved that the truly associated variables will be among those having the highest marginal coefficients.
In our simulation study, we only handle the unfaithfulness involving two associated variables X 1 and X 2 by using b 1 > 0, b 2 < 0 and r(X 1 , X 2 ) > 0 as illustrated in Figure 1(d). The marginal coefficientsβ 1 andβ 2 will be small due to the cancelation given by Equation (2). When b 1 > 0, b 2 > 0 and r(X 1 , X 2 ) < 0, the unfaithfulness also happens. This corresponds to a situation that the minor alleles of both X 1 and X 2 increase the diseases risk but X 1 and X 2 are negatively correlated, as illustrated in Figure 1(c). The simulation result shows that the marginal test and Lasso perform poorly. The better performance of BEAM should be attributed to its first order Markov chain designed for the accommodation of correlation. Although our exhaustive method performs reasonably well, the direct extension of our method to deal with three or more loci is computationally expensive. Therefore, solving the unfaithfulness issue is still challenging.

The Relationship between interaction models and unfaithfulness
In this work, we only deal with a two-locus association pattern involving the unfaithfulness. There are many works [3,5,6,30] discussing two-locus associations. Most of them belong to the category of interaction analysis (see details in [1,2]). The SNP interaction is also referred to as "epistasis". The most common statistical definition of interactions is the statistical deviation from the additive effects of two loci on the phenotype [2]. Using the same example we discussed in the introduction section, testing interactions between X 1 and X 2 is to first fit the regression model (or logistic regressions for case-control data) Y~b 1 X 1 + b 2 X 2 + b 12 X 1 X 2 and then test the significance of b 12 . There is no direct connection between b 1 (or b 2 ) and b 12 In the analysis of unfaithfulness, the relationship between marginal coefficients (β 1 ,β 2 ) and joint coefficients (b 1 , b 2 ) is given in Equation (2). The interaction term plays no role here. Therefore, it is not appropriate to use interaction models to detect associations masked by unfaithfulness.

The Relationship between unfaithfulness and confounding
Suppose we are studying the relationship between two variables X and Y using model Y ∼βX. Confounding arises when there is another observed variable Z which is independently associated with X and Y. Specifically, we haveβ xz = 0 for model X ∼β xz Z andβ yz = 0 fot model Y ∼β yz Z. When studying the relationship between X and Y, it is necessary to account for the confounding effect by using model Y~b yz Z + b yx X. In other words, confounding is more like the situation illustrated in Figure 1 (b). Readers are referred to [31] for the detailed explanation of confounding. The unfaithfulness is different. For model Y ∼β 1 X 1 and Y ∼β 2 X 2 , bothβ 1 andβ 2 are close to zero. For joint model Y~b 1 X 1 + b 2 X 2 , both b 1 and b 2 are not zero, as illustrated in Figure 1(c) and 1(d).

Biological interpretations
There are two possible biological interpretations. The first interpretation is illustrated in Figure 1(d). Consider two loci X p and X q which are positively correlated. When X q increases the disease risk (b q > 0) and X p acts as a protective locus (b p < 0), unfaithfulness happens. The identified associations and their coefficients listed in Table 3 indicate that these associations indeed exist.
The second interpretation is illustrated in Figure 1(c). Consider two loci X p and X q which are negatively correlated. When both X p and X q increase the disease risk (b p > 0 and b q > 0), unfaithfulness also happens. This case may be particularly interesting when analyzing SNPs with low allele frequencies [32]. Suppose the allele frequencies of both X p and X q are low and thus the mutations happening at these two loci are relatively recent. We can further assume the haplotype aa does not exist (because the probability of both two mutations happen in a short period is very small). This implies these two loci are negatively correlated. Unfortunately, we do not identify this type of associations. Possible reasons include: (1) The current genotyping chip is designed based on the "common disease/common variant" model [33,34], the low frequency SNPs are not directly assayed. (2) The statistical power of current testing strategy is relatively low to handle rare variants.

The unfaithfulness implications on tagSNPs
GWAS is considered as a powerful approach to detecting genetic susceptibility of common diseases. Such studies require the genotypes of a huge number of SNPs across the genome, each of which is tested for association with the phenotype of interest. This is generally referred to as the direct test of association, in which the functional mutation is presumed to be genotyped. An alternative approach is to take advantage of the correlation between SNPs. This approach genotypes a subset of SNPs, called tagSNPs, which are in high linkage disequilibrium with other SNPs [33]. The association tests of tagSNPs are used to indirectly infer the association of other correlated SNPs. This approach is widely used to save genotyping costs in GWAS. Many tagging methods [33,35,36] have been developed to reduce the number of markers to be genotyped. One key assumption in these methods is that the association analysis of a set of highly correlated SNPs is equivalent with the association analysis of tagSNPs of this set. However, the existence of unfaithfulness poses a challenge for these methods. The weak marginal association of a tagSNP does not imply the weak association of the corresponding genome region in which this tagSNP is located. The reason is that some non-genotyped SNPs correlating with the tagSNPs may jointly display strong associations in the presence of unfaithfulness.
In this work, we analyzed the WTCCC data generated by the Affymetrix 500 K chip and a Crohn's disease data set generated by the Illumina chip. The Affymetrix 500 K chip spaces SNPs along the genome and the Illumina chip uses the tagSNP design. LD becomes less apparent in the Illumina data set and we did not find any association masked by unfaithfulness. This result suggests that it is very difficult to detect these associations by using the tagSNP design. If more SNPs could be genotyped in the future GWAS, we would detect more unknown associations.

Conclusion
The phenomenon named "unfaithfulness" has been discussed as a mathematical concept in the literature. In this work, we answered the question whether there exist associations masked by unfaithfulness in genome-wide association studies. We developed a simple and fast method to examine all SNP pairs and demonstrated that our method is applicable to analyze genome-wide SNP data sets. We conducted experiments on both simulated data and seven real data sets from WTCCC and identify many associations masked by unfaithfulness. As expected, these identified associations only occur in local area. This implies that only the local search is needed to find such associations.
To date, we can only connect some identified associations to publicly available results from other association studies. As independent data set is limited as this moment, we have difficulties to replicate these findings. The biological interpretation of many findings remains unclear. It would be of great interest if their biological functions could be investigated. In addition, we only handle the two-locus associations in the presence of unfaithfulness. Detecting such associations for three or more loci is still an open problem.

Methods
Given a data set with L SNPs and n samples, we use X l to denote the l-th SNP, l = 1,···, L and Y to denote the class label (0 for control and 1 for case). SNPs are biallelic genetic markers. Capital letters (e.g. A, B,...) and lowercase letters (e.g. a, b,...) are often used to denote major and minor alleles, respectively. For simplicity, we use {0, 1, 2} to represent the three genotypes {AA, Aa, aa}, respectively.

Definition of the association masked by unfaithfulness
Considering a pair of loci X p and X q , four logistic regression models are typically involved to test associations masked by unfaithfulness: = β 0 +β p,1 I(X p = 0) +β p,2 I(X p = 1), (4) = β 0 +β q, 1 I(X q = 0) +β q, 2 I(X q = 1), (5) and M1⊕2 : log where I(V = v) is the indicator function that takes the value 1 if V = v is true and 0 otherwise. In order to make the representation of both logistic regression models and log-linear models (introduced later) in a compact and consistent way, we use the notation adopted in [37] and rewrite the above logistic regression models in the following forms: and Please note that the superscripts X p and X q in Equation (8), Equation (9) and Equation (10) are merely the labels and do not represent the exponents. The term β X p i represents the coefficient of X p at category i. Throughout this paper, we use ℳ to denote logistic regression models. We will use M to denote log-linear models. The log-likelihood function of a logistic model ℳ is denoted as L ℳ , and its maximum likelihood estimation (MLE) is denoted asL M . The log-likelihood function of a log-linear model M is denoted as L M , and its maximum likelihood estimation (MLE) is denoted asL M . For example, ℳ 1 is a logistic regression model whose log-likelihood function and MLE are denoted by L ℳ1 andL M 1 .
Our goal is to test if ℳ 1⊕2 is significantly different from ℳ 0 when both ℳ 1 and ℳ 2 are not. The likelihood ratio test is often used to conduct such tests. To test the difference between ℳ 1⊕2 and ℳ 0 , the following three steps are involved: 1. Fit a logistic regression model defined in Equation (10) and obtain the log-likelihoodL M 1 ⊕ 2 . 2. Compute the log-likelihoodL M 0 of the null logistic regression model defined in Equation (7). 3. Calculate P-value using the c 2 test on the value 2 Similarly, the test of difference between ℳ 1 (or ℳ 2 ) and ℳ 0 involves the following three steps: 1. Fit a logistic regression model defined in Equation (8) (or Equation (9)) to measure the main effect of X p (or X q ) and obtain the log-likelihoodL M 1 (or L M 2 ). 2. Compute the log-likelihoodL M 0 of the null logistic regression model defined in Equation (7). 3. Calculate P-value using the c 2 test on the value 2 Directly using regression methods for testing all pairs of SNPs in GWAS would be very time-consuming. Often the parallel computation was recommended [38]. Here, we propose to use log-linear models [37] instead of logistical regression models in GWAS. We show that this makes the likelihood ratio test computationally more efficient in genome-wide SNP data analysis. In the following, we briefly summarize the key components. The details are explained in the supplementary document (Additional file 1).

Likelihood ratio tests using log-linear models
Given two loci X p and X q , a contingency table of X p , X q , Y will be used to test the association masked by unfaithfulness between (X p , X q ) and Y. The size of the contingency table is I × J × K, where I = 3, J = 3, K = 2. We use n ijk to denote the observed count in the cell (i, j, k) in the contingency table (Table 4). Here n ijk is considered as the realization of a random variable N ijk assumed as Poisson-distributed in log-linear models.
We use the dot convention to indicate summation over a subscript, e.g., n i.. = j,k n ijk is the number of observations with X p = i. Similarly, we have n .j. = i,k n ijk and n ..k = i,j n ijk. We also have n ij. = k n ijk , n .jk = i n ijk and n i.k = j n ijk. Throughout this paper, we use μ M ijk to denote the expectation of N ijk under log-linear model M, and useμ M ijk to denote the MLE of μ M ijk . The equivalence between log-linear models and logistic models are given in Table 5 and the devianceL M 1 −L M 0 (orL M 2 −L M 0 ) can be computed aŝ In Equation (11) and Equation (12),μ P ijk andμ B ijk have the closed-form solutions (please check the supplementary document (Additional file 1) for the derivation): and = 0, due to multi-colli-  this solution works well in practice (We have compared our results with the standard software R, in which the multi-collinearity problem is elegantly handled when fitting generalized linear models. It turns out that our results agree with the results given by R). Then the test statistics can be efficiently computed. As a result, we are able to test every pair of loci to search for associations masked by unfaithfulness in GWAS. Table 6 gives the running time of our method for data sets of different sizes.
An exhaustive approach to detecting the two-locus associations masked by unfaithfulness in GWAS This approach involves the following steps: • Step 1. For all of L SNP markers, we first filter out those SNPs with significant main effects using Equation (12) since we are only interested in those markers without significant main effects. The L P-values can be adjusted by either the classic Benjamini-Hochberg method for controlling false discovery rate (FDR) or the Bonferroni correction for controlling family wise error rate (FWER).
• Step 2. For the remaining L' SNPs without significant main effects, we check every pair using the Equation (11). Again, the P-values can be adjusted for controlling either FDR or FWER.
The P-value thresholds in both Step 1 and Step 2 are specified by users. The default setting of the threshold is 0.1. The multiple factor for Bonferroni correction is L'(L' -1)/2, where L' is the number of SNPs after removing those SNPs with significant marginal effects.
Since the number of SNPs with significant marginal effects only accounts for a small fraction of the entire SNP set, we have L'(L' -1)/2 ≈ L(L -1)/2.

Simulation design
Let p(D|G i ) denote the probability of an individual being affected given its genotype G i (i.e., the penetrance of G i ), and let p(D|G i ) denote the probability of an individual not being affected given its genotype G i . Based on the definition of the odds of a disease the penetrance p(D|G i ) of the genotype G i can be calculated using The disease prevalence p(D) and genetic heritability h 2 are given as The odds table of our simulation model is given in Table 7. It is a multiplicative model of odds ratio, i.e., it is an additive model on the log-odds scale. The reason we choose this model is that we try to exclude interference of the interaction effect when we discuss the unfaithfulness. Essentially, the unfaithfulness arises due to the correlation cancelation. The interaction effects play no role here.
For simplicity, we restrict θ 11 = θ 12 = θ a and θ 21 = θ 22 = θ b . The parameter θ a > 1 means that the minor allele "a" increases the disease risk. This corresponds to the bivariate regression coefficient b 1 > 0. Similarly, θ b < 1 implies b 2 < 0. In the presence of linkage disequilibrium (linkage disequilibrium measure Δ > 0), unfaithfulness arises. To simulate this situation, we further set θ a = θ and θ b = 1/θ. In the simulation, the prevalence p(D) and the heritability h 2 are controlled by the parameters a and θ. We first specify the disease prevalence p(D) and genetic heritability h 2 . Then we numerically solve the parameters (a and θ) based on the above equations. We set p(D) = 0.1 and h 2 = 0.02. The details are given in the supplementary document (Additional file 1).

Additional material
Additional file 1: In the supplementary document (Additional le 1), we present the details of simulation. We also give a brief introduction to log-linear models which are used in the main article. Finally, we provide full lists of the results identified from the WTCCC data sets.