Abstract
Background
Genes in a functional pathway can have complex interactions. A gene might activate or suppress another gene, so it is of interest to test joint associations of gene pairs. To simultaneously detect the joint association between disease and two genes (or two chromosomal regions), we propose a new test with the use of genomic similarities. Our test is designed to detect epistasis in the absence of main effects, main effects in the absence of epistasis, or the presence of both main effects and epistasis.
Results
The simulation results show that our similarity test with the matching measure is more powerful than the Pearson's χ^{2 }test when the disease mutants were introduced at common haplotypes, but is less powerful when the disease mutants were introduced at rare haplotypes. Our similarity tests with the counting measures are more sensitive to marker informativity and linkage disequilibrium patterns, and thus are often inferior to the similarity test with the matching measure and the Pearson's χ^{2 }test.
Conclusions
In detecting joint associations between disease and gene pairs, our similarity test is a complementary method to the Pearson's χ^{2 }test.
Background
Genes in a functional pathway can have complex interactions. A gene might activate or suppress another gene, so it is of interest to test joint associations of gene pairs. Differing from epistasis (generally defined as the interaction between different genes [1]), joint associations herein include both main effects and interactions. Haplotypes from two receptors can trigger significant interactions affecting disease status [2]. Moreover, detecting associations with the use of haplotypes constructed by several adjacent and highly correlated singlenucleotide polymorphisms (SNPs) is an economical strategy. These all enlighten us regarding ways to develop methods for discovering gene pairs in association with disease by using haplotypes.
There is a growing interest in detecting genegene interactions [1,3,4], and some methods have been proposed to detect interactions. A wellknown approach to detecting SNPSNP interactions, the multifactor dimensionality reduction (MDR) method [58], however, has not been developed for testing haplotypehaplotype interactions. Another commonly used method is the classification and regression trees (CART) [912]. This concept has been extended to analyze haplotype data, known as the HapForest approach [13].
In this paper, we do not focus only on interactions because the definition of independence between two genes is arbitrary, often varying according to the field under discussion, such as biology, statistics or epidemiology [1]. Instead, we focus on detecting joint associations. To simultaneously detect joint association between disease and two genes (or two chromosomal regions), we propose a new test with the use of genomic similarities. Similaritybased methods are less vulnerable to the penalty of testing many markers or haplotypes, and can be more powerful than conventional association methods in some situations [14]. Our proposed test is designed to detect epistasis in the absence of main effects, main effects in the absence of epistasis, or the presence of both main effects and epistasis. We further compare our method with the HapForest approach [13], the Pearson's χ^{2 }test, and the tests for SNP × SNP epistasis via simulation studies.
Methods
Similarity Measures
Let
A. Diplotype perspective
A.1. Similarity measure based on identitybystate (IBS) allele sharing (referred to as 'IBS'):
where L is the number of loci considered in G_{k};
A.2. Similarity measure based on IBS inversely weighted by genotype frequencies (referred to as 'WIBS'):
where
Joint Similarity Regarding Two Genes
A similarity measure accounting for the joint association of genes G_{1 }and G_{2 }for the i^{th }and j^{th }subjects is
where
B. Haplotype perspective
B.1. Similarity based on the counting measure for haplotypes (referred to as 'COUNT'):
Let h_{i }and h_{j }be the i^{th }and j^{th }categories of haplotypes in a gene,
where
B.2. Similarity based on the matching measure for haplotypes (referred to as 'MATCH'):
Let h_{i }and h_{j }be the i^{th }and j^{th }categories of haplotypes in a gene, then the similarity based on the matching measure for haplotypes is
where s(h_{i}, h_{j}) is 1 only when all alleles match for the i^{th }and j^{th }haplotypes, otherwise s(h_{i}, h_{j}) is 0.
Joint Similarity Regarding Two Genes
Let
where
Similarity Test
Let the dissimilarity accounting for the joint association of genes G_{1 }and G_{2 }for the i^{th }and j^{th }subjects be
where n_{CS }and n_{CN }are the numbers of cases and controls, respectively;
Simulation Study
Simulation studies were conducted to evaluate the performance of our method. We extended the simulation scheme of Li et al. [18] to two chromosomal regions. In each region, 4,000 haplotypes across 300 kb were generated using the coalescentbased program ms [19]. The effective population size was set at 10,000, the recombination rate per base pair (bp) per generation was set at 10^{9}, and 300 SNPs were simulated in each region. For the human genome, recombination occurs at an average rate of about 10^{8 }per bp per generation [20]. Our recombination rate, 10^{9 }per bp per generation, is the low end of the recombination rates in the human genome [18], representing a stronger linkage disequilibrium (LD). We chose this rate because multimarker approaches are primarily designed for strongLD regions. In each chromosomal region, 2,000 diplotypes were generated by randomly pairing the 4,000 haplotypes. Then the 2,000 diplotypes of the first region were randomly paired with the 2,000 diplotypes of the second region, to form 2,000 subjects. In this way, we generated 300 datasets.
We then considered nine disease models listed in Additional file 1. Additional file 1 lists the causal allele frequencies, the penetrance values of twolocus genotypes, and the marginal penetrance values of onelocus genotypes, for all disease models. Model 0 was used to evaluate TypeI error rates, while the other eight models were used to evaluate powers. Models 16 exhibit interactions in the absence of main effects when genotypes conform to HardyWeinberg equilibrium. We used these six disease models because they further challenged the ability of our method to discover the joint associations (or 'interactions' in this situation) of gene pairs. Models 7 and 8 exhibit both interactions and main effects. Model 7 is the jointly dominantdominant model, which requires at least one copy of the disease allele from both loci to be affected [21,22]. Model 8 has the same penetrance table with Model 3, but has different causal allele frequencies. We deliberately let the causal allele frequency of one locus be smaller than that of another locus.
Additional file 1. Table S1. The penetrance tables and causal allele frequencies of nine disease models.
Format: DOC Size: 72KB Download file
This file can be viewed with: Microsoft Word Viewer
For each dataset, we first randomly selected two SNPs (each from among 300 SNPs in a region) with similar MAFs to those of the causal SNPs (the tolerable difference was set to be 0.02), pretending them as the two causal SNPs. We then used the Hclust method [23,24] to choose tag SNPs with a subset formed by 200 subjects randomly drawn from the pool of 2,000 subjects. Tag SNPs were chosen with quality (MAF > 0.1) and correlation (the cutoff value for finding clusters was set to be 0.85). In each repetition, cases and controls were sampled with replacement from the pool of 2,000 subjects, where case/control status was assigned according to the genotypes of the two causal SNPs. After generating the phenotypes, the genotypes of the causal SNPs were removed from our datasets. Each chromosomal region was formed by eight SNPs  four to the left and four to the right of every causal SNP.
We evaluated the performance of our method with the matching measure ('MATCH') and the counting measure ('COUNT') of haplotypes. We also used two genotype similarity measures: 'IBS' and 'WIBS'. We compared these with the HapForest approach [13]. HapForest is based on a tree structure, and is naturally suitable for analyzing interactions. Following the instructions of HapForest, we first invoked SNPHAP [25] to estimate the haplotype frequencies for each individual. Then HapForest was used to identify haplotypes and haplotypehaplotype interactions in association with the disease. This method suggests potential epistasis among significant haplotypes. For HapForest, a rejection of null hypothesis was defined as the identification of at least one significant haplotype from any of the two chromosomal regions.
The Pearson's χ^{2 }test was also performed for comparison, in which the joint haplotype distributions of the two chromosomal regions were compared between cases and controls. Rather than using the asymptotic χ^{2 }distribution, we randomly assigned the disease status in each permutation and determined the P value of observed χ^{2 }statistics. To calculate haplotype similarities from unphased multimarker genotypes, we first inferred haplotype phases by the EM algorithm, using the function of 'haplo.em' in the 'haplo.stats' package [17]. The obtained posteriors were then treated as weights, and all possible haplotype pairs were considered with their probabilities (see equation (2)). All the haplotypes with frequencies less than 0.01 are considered to be rare haplotypes. To avoid possible genotyping errors, we follow Sha et al. [26] to merge each rare haplotype with its most similar common haplotype (see the modified EM algorithm proposed by Sha et al. [26]). For example, Haplotype A (11121111) is considered to be a rare haplotype because its frequency is less than 0.01. Haplotypes C (11111111) and F (11122111) are the most similar haplotypes to Haplotype A (both with a similarity of 0.875 by using the counting measure), and their haplotype frequencies are 0.2 and 0.1, respectively. We merge Haplotype A with Haplotype C, the most similar haplotype with the highest frequency. We then update the haplotype data by replacing Haplotype A with Haplotype C.
We also compared our methods with the tests for SNP × SNP epistasis by using casecontrol
data or caseonly data (with the fastepistasis command implemented by PLINK1.07)
[27]), hereafter referred to as 'CSCN' and 'CS', respectively. In our simulation, each
chromosomal region was formed by eight SNPs, and there were 64 tests for SNP × SNP
epistasis. We recorded the minimum P value (P_{min}) from among all the 64 P values, and then adjusted this P_{min }on the basis of Sidak correction [28], with an effective number of tests, M_{eff}. That is, we adjusted the minimum P value (P_{min}) by
We then evaluated the validity and power of the eight tests with the 300 datasets. For each dataset, we recorded the P values of 50 repetitions (so there were 15,000 P values in total); in each repetition, P values were obtained with 1,000 permutations. Given a significance level, the type I error rate (if under Model 0) or power (if under Models 18) was the proportion of the number of P values smaller than the significance level to the total number of P values.
For CSCN and CS, the P value used was
Results
TypeI Error Rates
In Additional file 1, Model 0 (disease status independent of the composite genotypes) was used to evaluate the typeI error rates. This model demonstrates our null hypothesis: no main effects and no interactions. In this model, the penetrance of each composite genotype was set to be 0.05. The sample size was set at 200 subjects, of which half were cases and half were controls. Figure 1 presents the typeI error rates under different nominal significance levels (α). For α smaller than 0.2, the typeI error rates of all the tests corresponded to the nominal significance levels (α), suggesting the validity of these tests. (For α larger than 0.2, the typeI error rates of HapForest failed to match with the nominal significance levels. HapForest reported P values as 1.0 when the association signal was not strong. However, this makes no influence on our following discussions because α is usually set at a small value.)
Figure 1. TypeI error rates under different nominal significance levels. The xaxis is nominal significance level, and the yaxis is typeI error rate.
Statistical Power
For all models except for Models 2 and 7, the total sample size was set at 1,000 subjects, of which half were cases and half were controls. For Models 2 and 7, the total sample size was set at 150 and 50, respectively. If the sample size was also set at 1,000 for Models 2 and 7, the powers of these tests would be all close to 1. Therefore, we chose two smaller sample sizes for effectively exploring the power difference between these tests. The power performances of these tests vary with the property of disease mutants introduced at rare/common haplotypes.
We first define two scores to distinguish the two situations. Let
Figure 2 presents the powers of the eight tests when α is set to be smaller than 0.1, stratified by the property of disease mutants introduced at rare/common haplotypes. For most models, the two most powerful tests are our similarity method with the matching measure (MATCH) and the Pearson's χ^{2 }test. MATCH is more powerful than the Pearson's χ^{2 }test when the disease mutants were introduced at common haplotypes. Conversely, MATCH is less powerful than the Pearson's χ^{2 }test when the disease mutants were introduced at rare haplotypes. For Model 1, haplotypeperspective methods provide no power, while diplotypeperspective methods (IBS and WIBS) and the test for SNP × SNP epistasis by using caseonly data (CS) have better performances.
Figure 2. Powers of the eight tests, stratified by the property of disease mutants introduced at rare/common haplotypes. The xaxis is significance level, and the yaxis is power. The top row is for disease mutants introduced at rare haplotypes; the bottom row, at common haplotypes. The numbers shown in the parentheses are the numbers of repetitions summed from all the datasets with disease mutants introduced at rare/common haplotypes.
HapForest is not as powerful as MATCH and the Pearson's χ^{2 }test. HapForest suggests potential epistasis among significant haplotypes. At each step, it builds a classifier that optimally distinguishes cases from controls based on haplotype data. This divides the whole sample into smaller and smaller subgroups by maximizing the local optimality at each node. However, the combination of local optimalities does not assure us of an overall optimality [32].
The tests for SNP × SNP epistasis by using casecontrol data or caseonly data (CSCN and CS) are not powerful under most disease models. Although our disease status was influenced by the joint effects of two SNPs (see Additional file 1), the tests for SNP × SNP epistasis suffered from power loss because of the need of corrections for multiple testing.
COUNT and IBS often have similar performances, because similarity measure based on the number of alleles in common between haplotypes (COUNT) is similar to that based on the number of alleles in common between individuals (IBS). Model 1 is an exception, because haplotypeperspective methods would not present any power under this model (see the penetrance values of twolocus genotypes for Model 1). WIBS is a counting measure inversely weighted by genotype frequencies, and it is more powerful than COUNT and IBS. For most models (Models 26 and 8), COUNT, IBS, and WIBS are inferior to MATCH, because the counting measures are more sensitive to marker informativity and LD patterns (results not shown). For Model 7, it requires at least one copy of the disease allele from both loci to be affected. Because the disease status is influenced by the counts of disease alleles, methods with the counting measures (COUNT, IBS, and WIBS) are more powerful.
Discussion
Detecting joint associations of candidate genes responsible for common human diseases is a wellrecognized issue. A candidate gene can contain many SNPs, and highdimensionality becomes an important issue. The Pearson's χ^{2 }test and the tests for SNP × SNP epistasis suffer from power loss because of large numbers of degrees of freedom and the need of adjustment for multiple testing, respectively. Compared with these conventional association methods, similarity methods are less vulnerable to the penalty of highdimensionality.
Some similarity methods have been proposed based on this consideration. Tzeng et al. [15] compared the casecase similarity with the controlcontrol similarity, because haplotypes around a causal locus might be more similar in two cases than in two controls randomly selected from the population. However, as pointed out by Sha et al. [26], this consideration might not be very plausible for complex diseases which were presumed to be affected by many genes and geneenvironment interactions. The similarity within controls is not necessarily smaller than that within cases, because controls could be more likely to share protective haplotypes. Therefore, Sha et al. [26] proposed a test statistic that compared the betweengroup similarity with the withingroup similarity. Our test statistics is also based on this consideration. Our test and Sha et al.'s test [26] will have similar performances, given a same similarity measure.
In this paper, we use the product of similarities of two genes/regions as a new similarity
measure, which can account for the joint association of the two genes/regions, including
main effects and/or interactions. This new measure can be built in the similarity
test statistic. Furthermore, our equation (3) can be used to test the main effects
(see Appendix II) or the joint associations of gene triplets by using:
The computational burden of our method is reasonable for real data analyses, although permutation is required to obtain P values. If there are 100 candidate genes (each with eight tag SNPs), there will be a total of 4,950 combinations of gene pairs. With our experiences in simulations, it might take two to three days to test the 4,950 combinations for approximately 1000 subjects, given an Intel Xeon workstation with four 2.0 GHz CPUs and 2.0 GB of memory.
In general, our similarity test with the matching measure (MATCH) and the Pearson's χ^{2 }test have better power performances. However, because both are haplotypeperspective methods, they are not appropriate for Model 1. Under this model, only the four heterozygous genotypes (AABb, AaBB, Aabb, aaBb) lead to the disease. The implication is that besides the withinlocus interference, there is some betweenlocus interference, and the two interferences cancel out [21] (so the doubleheterozygosity genotype does not lead to the disease). The four heterozygous genotypes (AABb, AaBB, Aabb, aaBb) generate four combinations of haplotypes: AB (one with allele A and one with allele B), Ab, aB, ab, with a same probability. Therefore, the four combinations of haplotypes are equally distributed in cases and in controls, and the haplotypeperspective methods cannot provide any power to this model.
The concept of testing joint associations can be used in the genomic distancebased
regression [16]. Let D be the distance/dissimilarity matrix with elements:
Conclusions
In detecting joint associations between disease and gene pairs, our similarity test is a complementary method to the Pearson's χ^{2 }test.
Competing interests
The authors declare that they have no competing interests.
Authors' contributions
WY L conceptualized the study, performed the simulation studies, and drafted the manuscript. WC L provided advice and revised the manuscript. All authors read and approved the final manuscript.
Appendix
Appendix I: Derivation of equation (3)
Similarly,
Note that
Appendix II: Test for main effects
When testing for main effects, we use only one gene/region in equation (3), i.e.,
where
Acknowledgements
We thank the three anonymous reviewers for their constructive comments. This study was partly supported by National Science Councils, Taiwan.
References

Cordell HJ: Epistasis: what it means, what it doesn't mean, and statistical methods to detect it in humans.
Hum Mol Genet 2002, 11(20):24632468. PubMed Abstract  Publisher Full Text

Lin M, Li H, Hou W, Johnson JA, Wu R: Modeling sequencesequence interactions for drug response.
Bioinformatics 2007, 23(10):12511257. PubMed Abstract  Publisher Full Text

Cordell HJ: Estimation and testing of geneenvironment interactions in familybased association studies.
Genomics 2009, 93(1):59. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Musani SK, Shriner D, Liu N, Feng R, Coffey CS, Yi N, Tiwari HK, Allison DB: Detection of gene × gene interactions in genomewide association studies of human population data.
Hum Hered 2007, 63(2):6784. PubMed Abstract  Publisher Full Text

Hahn LW, Ritchie MD, Moore JH: Multifactor dimensionality reduction software for detecting genegene and geneenvironment interactions.
Bioinformatics 2003, 19(3):376382. PubMed Abstract  Publisher Full Text

Moore JH, Gilbert JC, Tsai CT, Chiang FT, Holden T, Barney N, White BC: A flexible computational framework for detecting, characterizing, and interpreting statistical patterns of epistasis in genetic studies of human disease susceptibility.
J Theor Biol 2006, 241(2):252261. PubMed Abstract  Publisher Full Text

Ritchie MD, Hahn LW, Moore JH: Power of multifactor dimensionality reduction for detecting genegene interactions in the presence of genotyping error, missing data, phenocopy, and genetic heterogeneity.
Genet Epidemiol 2003, 24(2):150157. PubMed Abstract  Publisher Full Text

Ritchie MD, Hahn LW, Roodi N, Bailey LR, Dupont WD, Parl FF, Moore JH: Multifactordimensionality reduction reveals highorder interactions among estrogenmetabolism genes in sporadic breast cancer.
Am J Hum Genet 2001, 69(1):138147. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Breiman L: Classification and regression trees. Belmont, CA: Wadsworth International Group; 1984.

Bureau A, Dupuis J, Falls K, Lunetta KL, Hayward B, Keith TP, Van Eerdewegh P: Identifying SNPs predictive of phenotype using random forests.
Genet Epidemiol 2005, 28(2):171182. PubMed Abstract  Publisher Full Text

Clark LA, Pregibon D: Treebased models. In Statistical Models in S. Edited by Chambers JM, Hastie TJ. Pacific Grove, California: Wadsworth and Brooks/Cole Advanced Books and Software; 1992:377419.

Zhang H, Bonney G: Use of classification trees for association studies.
Genet Epidemiol 2000, 19(4):323332. PubMed Abstract  Publisher Full Text

Chen X, Liu CT, Zhang M, Zhang H: A forestbased approach to identifying gene and gene gene interactions.
Proc Natl Acad Sci USA 2007, 104(49):1919919203. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Lin WY, Schaid DJ: Power comparisons between similaritybased multilocus association methods, logistic regression, and score tests for haplotypes.
Genet Epidemiol 2009, 33(3):183197. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tzeng JY, Devlin B, Wasserman L, Roeder K: On the identification of disease mutations by the analysis of haplotype similarity and goodness of fit.
Am J Hum Genet 2003, 72(4):891902. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wessel J, Schork NJ: Generalized genomic distancebased regression methodology for multilocus association analysis.
Am J Hum Genet 2006, 79(5):792806. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Schaid DJ, Rowland CM, Tines DE, Jacobson RM, Poland GA: Score tests for association between traits and haplotypes when linkage phase is ambiguous.
Am J Hum Genet 2002, 70(2):425434. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Li Y, Sung WK, Liu JJ: Association mapping via regularized regression analysis of singlenucleotidepolymorphism haplotypes in variablesized sliding windows.
Am J Hum Genet 2007, 80(4):705715. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hudson RR: Generating samples under a WrightFisher neutral model of genetic variation.
Bioinformatics 2002, 18(2):337338. PubMed Abstract  Publisher Full Text

A haplotype map of the human genome
Nature 2005, 437(7063):12991320. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Li W, Reich J: A complete enumeration and classification of twolocus disease models.
Hum Hered 2000, 50(6):334349. PubMed Abstract  Publisher Full Text

Neuman RJ, Rice JP: Twolocus models of disease.
Genet Epidemiol 1992, 9(5):347365. PubMed Abstract  Publisher Full Text

Rinaldo A, Bacanu SA, Devlin B, Sonpar V, Wasserman L, Roeder K: Characterization of multilocus linkage disequilibrium.
Genet Epidemiol 2005, 28(3):193206. PubMed Abstract  Publisher Full Text

Roeder K, Bacanu SA, Sonpar V, Zhang X, Devlin B: Analysis of singlelocus tests to detect gene/disease associations.
Genet Epidemiol 2005, 28(3):207219. PubMed Abstract  Publisher Full Text

Clayton D: SNPHAP  A program for estimating frequencies of large haplotypes of SNPs. Department of Medical Genetics, Cambridge Institute for Medical Research, Cambridge; 2006.

Sha Q, Chen HS, Zhang S: A new association test using haplotype similarity.
Genet Epidemiol 2007, 31(6):577593. PubMed Abstract  Publisher Full Text

Purcell S, Neale B, ToddBrown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, et al.: PLINK: a tool set for wholegenome association and populationbased linkage analyses.
Am J Hum Genet 2007, 81(3):559575. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Sidak Z: Rectangular confidence regions for the means of multivariate normal distributions.
Journal of the American Statistical Association 1967, 62:626633. Publisher Full Text

Cheverud JM: A simple correction for multiple comparisons in interval mapping genome scans.
Heredity 2001, 87(Pt 1):5258. PubMed Abstract  Publisher Full Text

Nyholt DR: A simple correction for multiple testing for singlenucleotide polymorphisms in linkage disequilibrium with each other.
Am J Hum Genet 2004, 74(4):765769. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Galwey NW: A new measure of the effective number of tests, a practical tool for comparing families of nonindependent significance tests.
Genet Epidemiol 2009, 33(7):559568. PubMed Abstract  Publisher Full Text

Briollais L, Wang Y, Rajendram I, Onay V, Shi E, Knight J, Ozcelik H: Methodological issues in detecting genegene interactions in breast cancer susceptibility: a populationbased study in Ontario.
BMC Med 2007, 5:22. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text