Abstract
Background
Although singleSNP analysis has proven to be useful in identifying many diseaseassociated loci, regionbased analysis has several advantages. Empirically, it has been shown that regionbased genotype and haplotype approaches may possess much higher power than singleSNP statistical tests. Both high quality haplotypes and genotypes may be available for analysis given the development of next generation sequencing technologies and haplotype assembly algorithms.
Results
As generally it is unknown whether genotypes or haplotypes are more relevant for identifying an association, we propose to use both of them with the purpose of preserving high power under both genotype and haplotype disease scenarios. We suggest two approaches for a combined association test and investigate the performance of these two approaches based on a theoretical model, population genetics simulations and analysis of a real data set.
Conclusions
Based on a theoretical model, population genetics simulations and analysis of a central corneal thickness (CCT) Genome Wide Association Study (GWAS) data set we have shown that combined genotype and haplotype approach has a high potential utility for applications in association studies.
Keywords:
Genotypebased tests; Haplotypebased tests; Association analysis; Test statistic combinationBackground
The development of genotyping and sequencing technologies has enabled scientists to investigate the impact of genomic loci on complex disorders and traits. Indeed, genomewide association studies (GWAS) and sequencing studies have identified many common singlenucleotide polymorphisms (SNPs) (for GWAS publication list, see http://www.genome.gov/gwastudies/ webcite) and rare variations [14] associated with common diseases. Although singleSNP analysis has proven to be useful in discovering many diseaseassociated loci, this strategy may be limited due to very stringent significance threshold and poor reproducibility [5]. Regionbased association studies have the advantages of less stringent significance level and potentially higher power if multiple associated variants are found within a region. Indeed, several empirical studies have demonstrated the superiority of genotype genebased association analysis over singleSNP strategy [6,7]. Also, there is some theoretical [8,9] and empirical evidence that haplotypebased tests may possess higher power than SNPbased tests. When intending to use haplotypes in an association study, one faces a problem of phase inference. While several statistical algorithms have been developed to infer unknown haplotypes from genotype data [1012], the improvements of sequencing technologies will enable researchers to assemble haplotypes from sequencing data with very high accuracy (for examples of existing assembly algorithms, see Bansal et al. [13], Bansal et al. [14], and Schatz et al. [15]). This opens up the opportunity to use highquality haplotypes and genotypes in sequencing association studies.
Numerous studies have reported cases when haplotypebased analysis resulted in detection of an association, while SNPbased analysis either did not yield any significant results or yielded much higher pvalues [1619]. A haplotypebased test may be more powerful than a genotypebased test if haplotypes tag a true causal variant better (although the imputation of untyped SNPs using publicly available reference panels may also be a powerful strategy), or if a SNPSNP interaction is present within a region. In general, it is unknown whether haplotype or genotypebased tests are more relevant for identifying an association of a genomic region with a phenotype. In this article we propose two statistical tests that explicitly combine both genotype and haplotype information for the purpose of preserving high power under both genotype and haplotype disease scenarios. We investigate two methods based on a combination of pvalues from genotype and haplotypebased association tests. The first method is a minimum of pvalues (MinPval), and the other is a sum test statistic based on inverse standard normal transformation of two pvalues (SumPval). Based on simulations, theoretical power calculations and application to a GWAS data set, we have highlighted the merits and the drawbacks of genotype and haplotypebased tests, and those of our combined approaches. The major conclusions from our work are as follows:
1. Combination of haplotype and genotypebased test statistics preserves power for both genotype and haplotype disease models;
2. In some of the considered scenarios, the performance of the MinPval approach is comparable to those of the SumPval method;
3. MinPval is much more robust than SumPval when one of the underlying tests has low power.
Methods
Genotype and haplotypebased tests
Let us assume that we are interested in testing the joint association of all the variants within a genomic region with either a dichotomous phenotype or quantitative trait. Next, assume we have chosen the two statistical tests for a regionbased association analysis: one genotype and one haplotypebased test. For haplotypebased tests haplotypes can be inferred from genotypes [1012] or assembled from sequencing data [13,15,20]. Several conventional genotypebased methods [2123] are applicable for common variants testing, whereas for sequencing data numerous recentlydeveloped rare variants approaches are available [2428]. Haplotypebased methodologies have also been extensively published elsewhere [2931], including rare haplotype tests [3234].
The combined approaches
Let us denote pvalues from a genotype and a haplotypebased tests as p_{1} and p_{2} respectively. Our first approach is SumPval [35]. Let us consider the inverse standard normal transformation of both pvalues and which are distributed as standard normal random variables under the null hypothesis. Here, we assume that y_{1},y_{2} is bivariate normal. The SumPval test statistic is P_{sum} = y_{1} + y_{2}. Under the null hypothesis, it is distributed as a normal random variable with zero mean and variance Var(y_{1} − y_{2}) = Var(y_{1}) + 2Cor(y_{1}, y_{2}) − Var(y_{2}) = 2 + 2p, where p is a correlation coefficient between y_{1} and y_{2}, since two statistical tests for the same genomic region may not be independent. The correlation coefficient p may be estimated via permutation procedure. The rejection region is large values of the test statistic, which is equivalent to low values for p_{1} and/or p_{2}. The theoretical pvalue for SumPval test is calculated as , where Φ(a, b, c) is a value of normal cumulative distribution function with mean b and variance c taken at the point a.
Our second approach (MinPval) is to utilize the minimum of the two pvalues as a test statistic, namely, min{p_{1}, p_{2}} [36]. Let us represent a test statistic as min{p_{1}, p_{2}} = min{1 − Φ(y_{1}), 1 − Φ(y_{2})}, where y_{1} and y_{2} defined above are distributed as standard normal random variables under the null. Thus, the theoretical cumulative distribution function of MinPval test statistic under the null hypothesis can be calculated as follows:
where 0 < x < 1. Given the rejection region is small values of the test statistic, theoretical pvalue for MinPval test is straightforward to compute using (1).
Theoretical power model
Within our theoretical framework the following model is adopted: the two test statistics S_{g} and S_{h} of the underlying genotype and haplotypebased tests, respectively, are assumed to asymptotically follow central chisquared distribution with 1 degree of freedom under the null hypothesis, and noncentral chisquared and with NCPs a and b, respectively, under the alternative hypothesis. One of the examples of the test which results in such null and alternative distributions is Rao’s score test on, for example, genotype or haplotype scores described in Additional file 1. Since the two tests are applied to the same data, the chisquared test statistics are likely to be positively correlated. The correlation between the two test statistics may vary from very low to high. For example, if within a region there are few SNPs in very high LD then we would expect the correlation between the tests to be high. Alternatively, we would expect the correlation to be low when variants within a region are independent. The correlation is modeled via underlying multivariate normal distribution, namely, to simulate the test statistics and a bivariate normal random vector with mean , unit variances and correlation coefficient p > 0 is generated, and the squares of the coordinates are taken as the proxy for the test statistics: , . To estimate the power of MinPval and SumPval tests we simulated 500,000 independent pairs (S_{g}, S_{h}) under the alternative huypothesis, calculated the test statistics for the combined approaches, and noted the share of statistically significant pairs. This procedure was done for every theoretical scenario (see “Results” section).
Additional file 1: Table S1. ShapiroWilk bivariate normality test pvalues for the genomewide significant genes, additional simulation and real data analysis results. ShapiroWilk test was used in real data analysis to justify our assumption of bivariate normality for calculation of theoretical pvalues for MinPval and SumPval tests. Additional simulations and real data analysis were performed using different pair of underlying tests.
Format: DOCX Size: 40KB Download file
Population genetics simulations
King et al. [37] provided the SFS_CODE (http://sfscode.sourceforge.net webcite) implementation of population genetics simulation for ANGPTL4 gene exons (http://home.uchicago.edu/~crk8e webcite). The authors assumed the demographic and distribution fitness effect parameters from Boyko et al. [38] and Gutenkunst et al. [39]. Using the SFS_CODE program 1000 haplotype pools each containing 20000 sampled “individuals” (40000 chromosomes) from a European population were generated. A data replicate was created from each haplotype pool by iterative random sampling of two haplotypes (thus, defining the genotype of an “individual”) and assigning a dichotomous phenotype conditional on a multisite genotype or a pair of haplotypes. Each data replicate contains 500 cases and 500 controls. Let us assume that there are L variants within the genomic region of interest, and the genotype of “an individual” is constructed from the sampled haplotypes. To describe the genotypebased disease model let us, without loss of generality, denote the genotypes at rare (MAF < 1%) causal SNPs as , causal common SNP g_{c1} (if present depending on a model), and other SNPs . Let us also define the assigned odds ratios of causal variants {b_{1}, …, b_{c − 1}}. The probability of a disease P(A) for a genotypebased scenario is calculated from the following:
For the haplotypebased scenarios let us consider the two sampled haplotypes {h_{1},h_{2}}. Also, denote H_{r} and H_{c} as the sets of rare (frequency in a haplotype pool <1%) and common causal haplotypes, respectively (depending on the disease model H_{c} may be empty). The probability of a disease is defined as:
where l{A} is an indicator of an event A, and {d_{r},d_{c}} are the odds ratios for causal rare and common haplotypes, respectively. For our simulations, we considered three phenotype models: “Rare” (only rare variants or haplotypes are riskcontributing), “Common” (only common variants or haplotypes), and “Both” (both types of variants or haplotypes). Following the scenarios of exomescale simulations of Wu et al. [40], we have assigned 50%, 20% and 10% of the observed rare variants (haplotypes) to be causal. Additionally, we chose one common causal SNP (haplotype) for “Both” and “Common” models. The odds ratios in (2) and {d_{r},d_{c}} in (3) were assigned as follows:
– for the “Rare” model: b_{l} = d_{r} = 4,l = 1,..,c and b_{c1} = d_{c} = 0
– for the “Both” model: b_{l} = d_{r} = 3,l = 1,..,c and b_{c1} = d_{c} = 1.2
– for the “Common” model: b_{l} = d_{r} = 1.5,l = 1,..,c and b_{c1} = d_{c} = 2
The average number of variants across data replicates is shown in Table 1.
Table 1. The average number of variants within a region across 1000 data replicates in population genetics simulations
Real data analysis
For the purpose of demonstrating the performance of the described methodologies we conducted a genebased analysis of the central corneal thickness (CCT) GWAS data sets described in Vithana et al. [41]. Briefly, the Singapore Indian Eye Study (SINDI), which is part of the Singapore Indian Chinese Cohort Eye Study (SICC) [42], consists of 2538 Indian subjects aged 40 and above, and the Singapore Malay Eye Study (SiMES) [4345] is a genomewide association study of CCT phenotype which contains 2542 Malay subjects aged 40 and above. Both SiMES and SICC adhered to the Declaration of Helsinki. Ethics approval for the both studies was obtained from the Singapore Eye Research Institute Institutional Review Board [41]. The combined data set consists of 5080 individuals genotyped at 552318 SNPs after quality control. In total, 5049 individuals were analyzed after excluding those with missing phenotype. Also, we attempted to replicate all the genomewide significant regions using Chinese samples from the SICC. This data set contains 2837 samples with nonmissing phenotype and covariates (age and gender). SNPs were mapped to genes based on the method outlined by Zhao et al. [46]. Briefly, information on gene identifiers (IDs), names, start and end positions on a chromosome were downloaded from the NCBI Genome database (http://www.ncbi.nlm.nih.gov/Genomes webcite). Gene regions included 10 kb upstream and downstream. Hierarchical mapping scheme (coding > intronic > 5′UTR > 3′UTR) was used if a variant was within 10 kb of multiple genes. The remaining intergene variants between two genes were grouped together. Haplotype inference was performed using the software Beagle [10] with reference panel from 1000Genomes Project (http://www.1000genomes.org/ webcite). In our analysis we adjusted for age, gender and the first ten principal components from Eigenstrat [47].
Statistical tests for population genetics simulations and real data application
Sequence Kernel Association Test (SKAT), introduced by Wu et al. [40], is a variance component score test derived from a semiparametric regression model. It was initially proposed to test the association of phenotype with multisite genotype; however, we also used SKAT to test an association of phenotype with haplotypes as described below. To show the consistency of empirical results we applied the same pair of underlying tests, namely, genotype SKAT and haplotype SKAT with linear kernel and uniform weights, to both population genetics simulations and real data. For genotype SKAT all rare variants (MAF < 1% in the sample) within a region were collapsed according to the method described by Thalamuthu et al. [48]. Briefly, the collapsed supervariant is the sum of minor alleles across rare variants within a region; if this sum is greater than 2 then the value of 2 is assigned. We did not apply weighting as described by Wu et al. [40] because that would substantially decrease power to identify an association with common SNPs. Haplotypebased SKAT is SKAT applied to a haplotype regression matrix R, which is constructed similar to those used by Zaykin et al. [31]. First, we pooled all rare haplotypes into one haplotype group, whereas each of the common haplotypes formed a separate haplotype group. Let us define the following notations: n is a number of individuals; {H_{1},…,H_{M}} is the haplotype groups with H_{M} being the most common group; is a haplotype regression matrix; and {h_{i,1};h_{i,2}} is a pair of haplotypes for i th individual. The haplotype matrix {R_{ij}} is constructed as follows:
where l{A} is an indicator of an event A. If there were no common haplotypes within a region, we formed three groups of haplotype: those with a frequency less that 0.05%, those in between 0.05% and 0.1%, and those with a frequency greater than 0.1%. For both tests we used the R (http://www.rproject.org/ webcite) package SKAT (http://www.hsph.harvard.edu/~xlin/software.html webcite). For population genetics simulations pvalues for all the tests were estimated using 1000 permutations. In real data analysis for the underlying tests we used theoretical pvalues as we believe that they reasonably approximate empirical pvalues given large sample size and normallydistributed quantitative trait [41]. Then we tested an assumption of bivariate normality by applying the ShapiroWilk test (R package “mvnormtest” http://cran.rproject.org/web/packages/mvnormtest/ webcite). If the normality test was not significant on the genomewide level, we used theoretical pvalues for both SumPval and MinPval; otherwise we used permutations. The permutation procedure and estimation of correlation coefficient are described in the next section.
Permutation procedure and estimation of correlation coefficient
To calculate theoretical pvalues for the proposed methods we estimated the correlation coefficient p using 500 permutations. The difficulty in applying permutations lies in the fact that the permutation procedure should preserve the relationship between all the covariates, and also between phenotype and covariates, but disrupt the relationship between phenotype and genotype. Several techniques have been developed for conducting permutation tests of partial coefficients in a multiple regression model [4951]. Among them the permutation of residuals under the reduced model [49] was shown to preserve correct type1 error for ttest [52] and was previously applied to microarray data analysis [53]. As the SKAT test can be obtained from a semiparametric regression model [40], let us consider the following genotype and haplotype regression models: Y = a_{1} − f_{1}(P) − Cc − ϵ and Y = a_{2} − f_{2}(R) − Cc − ϵ, where P is n × L collapsed genotype matrix, n is the sample size, L is the number of common SNPs within a region plus one for collapsed rare variants superlocus, Y is n × 1 vector of quantitative phenotype (CCT), C is n × 12 matrix of covariates which include age, gender and the first ten genotype principal components obtained from Eigenstrat [47], R is haplotype regression matrix, and f_{1} and f_{2} are unknown functions. To obtain the permutation values for the test statistics the reduced model Y = a − Cc − ϵ is fitted, and a, c, ϵ are the estimated constant coefficient, regression coefficients and residuals, respectively. Next, the residuals ϵ are permuted to obtain ϵ, and Y = a + Cc − ϵ. The permuted statistic values for both genotype and haplotype SKAT tests are calculated as respective SKAT statistics from semiparametric models Y = a_{1} − f_{1}(P) − Cc − ϵ and Y = a_{2} − f_{2}(R) − Cc − ϵ. Each pvalue obtained from permutations was transformed using the inverse standard normal transformation, and the value of p was estimated by a Pearson correlation coefficient.
Results
Theoretical power results
Depending on the disease model, one of the underlying tests (genotype or haplotypebased) is expected to be more powerful than the other underlying test. So, we assume that under the alternative hypothesis the noncentrality parameter (NCP) of the more powerful underlying test is , and the NCP of the less powerful underlying test is b = a/2. Figures 1 and 2 (Panel 1) show the power of MinPval and SumPval strategies as a function of correlation coefficient p and NCP a at the fixed type1 error of 0.05. As can be seen in Panel 1, power in general decreases slightly with increasing correlation. Panel 2 depicts the difference in power between the combined approaches and the more powerful underlying test. It is notable that both MinPval and SumPval achieved greater power than the more powerful underlying test for lower correlation. Also, MinPval approach lost a maximum of 5% power for high correlation and gained a maximum of 2% for low correlation, whereas for SumPval these values are more than 5% in both cases.
Figure 1. Performance of MinPval approach under the theoretical models. Different lines across the panels correspond to different levels of correlation p  0, 0.3, 0.6 and 0.9. Panel 1: Power of MinPval test as a function of NCP a and correlation p; Panel 2: Difference in power between MinPval and the more powerful underlying test as a function of NCP a and correlation p; Panel 3: Difference in power between MinPval and the less powerful underlying test as a function of NCP a and correlation p; Panel 4: Power of MinPval test as a function of correlation p and the ratio of NCP b to NCP of the more powerful underlying test (the latter is fixed at 10.5).
Figure 2. The performance of SumPval approach under the theoretical models. Different lines across the panels correspond to different levels of correlation p  0, 0.3, 0.6 and 0.9. Panel 1: Power of SumPval test as a function of NCP a and correlation p; Panel 2: Difference in power between SumPval and the more powerful underlying test as a function of NCP a and correlation p; Panel 3: Difference in power between SumPval and the less powerful underlying test as a function of NCP a and correlation p; Panel 4: Power of SumPval test as a function of correlation p and the ratio of NCP b to NCP of the more powerful underlying test (the latter is fixed at 10.5).
In Panel 3, where the difference in power between the combined approaches and the less powerful underlying test is shown, it can be seen that both MinPval and SumPval are consistently better than the less powerful test. This suggests that combination of statistical tests may prove beneficial when the underlying disease model is unknown. To investigate the impact of change of NCP b on the performance of the proposed approaches we fixed NCP a to be equal to 10.5 (corresponding to 90% power of a chisquared test with distribution under the alternative hypothesis, the type1 error is 0.05). Panel 4 of Figures 1 and 2 depicts the power of MinPval and SumPval as a function of correlation and a “fraction of NCP” – the ratio of b to 10.5. As can be seen in Panel 4, MinPval test achieved higher power than SumPval in the majority of scenarios. It is notable that SumPval lost much power when the value of b is low. Hence, MinPval approach is more robust with respect to underperformance of one of the underlying tests.
Population genetics simulation results
Panel 4 of Figure 3 shows the empirical type1 error estimate for the theoretical level of 0.05 for all the tests. The estimate of type1 error is distributed as a binomial random variable with 1000 trials and the probability of success 0.05 under the hypothesis of no inflation. The onesided 99% quantile of the described distribution is 0.067. As can be seen, in our simulations the type1 error was well controlled for all the tests.
Figure 3. Power comparison of genotypebased SKAT, haplotypebased SKAT, MinPval and SumPval tests for population genetics simulations, and an estimate of empirical type1 error. In each panel the top three disease models correspond to the haplotypebased disease scenario, whereas the lower three correspond to the genotypebased scenario. Disease models “Rare”, “Both” and “Common” are described in the section “Population genetics simulation”. Type1 error is set to 5%. Panel 1: 50% of rare variants/haplotypes were assumed to be causal; Panel 2: 20% of rare variants/haplotypes were assumed to be causal; Panel 3: 10% of rare variants/haplotypes were assumed to be causal; Panel 4: empirical type1 error estimate for simulations under the null hypothesis.
Panles 1–3 of Figure 3 depict the results of population genetics simulations analysis for all the phenotype models with 50%, 20% and 10% or rare causal variants/haplotypes, respectively, at the fixed 5% type1 error. For all the tests 1000 permutations were performed to estimate pvalues. Haplotypes were assumed to be known without ambiguity. Under the genotypebased disease scenarios, genotype SKAT is expected to be more powerful than haplotype SKAT, and vice versa under the haplotypebased scenarios. However, genotype SKAT was less powerful for many genotypebased phenotype models. A possible explanation of this observation is that when rare variants are strongly associated with phenotype, for some statistical tests pooling of rare haplotypes may be a better strategy than pooling of rare variants. Also, it should be noted that with the decrease in the percentage of causal rare variants/haplotypes, the power for “Rare” and “Both” phenotype models decrease substantially since for these models rare variants/haplotypes are the major carriers of an association signal. For “Common” phenotype model one common variant/haplotype has a significant impact on phenotype; so, the decrease in power with the lower proportion of causal rare variants/haplotypes is not as high as for other phenotype models.
As can be seen from the Panels 1–3 of Figure 3, for all the phenotype disease models, when both underlying tests were almost equally powerful (e.g. Panel 1 haplotype disease scenario “Common” model, and genotype disease scenario “Both” and “Common” models), the power of both MinPval and SumPval were on the same level or even higher than those of the underlying tests. However, when genotypebased SKAT significantly underperformed haplotypebased SKAT (e.g. Panel 1 haplotype disease scenario “Rare” and “Both”models), MinPval approach showed slightly lower power than the more powerful underlying test and greater power compared with SumPval approach. The maximum power loss of SumPval and MInPval compared with the more powerful underlying test across all phenotype models was 6.3% and 3.8% respectively (haplotype disease scenario “Both” model). These results are consistent with those obtained from the theoretical power considerations, and illustrate the great potential of the proposed methods in their application to real association studies.
To examine the effect of phasing on our results we repeated the analysis using the most probable haplotypes inferred by Beagle [10]. The reference panel consisted of 1094 simulated individuals to mimic the size of the publicly available reference panel from the 1000 Genomes Project (http://www.1000genomes.org webcite). The results of this analysis were very similar to those described above (data not shown). In addition, we applied the proposed methods with a different pair of underlying tests. The results are similar to those described above. For more details, see Additional files 1 and 2.
Additional file 2. Power comparison of the gene score haplotype test, the gene score genotype test, MinPval and SumPval statistical tests for population genetics simulations, and an estimate of empirical type1 error. In each panel the top three disease models correspond to the haplotypebased disease scenario, whereas the lower three correspond to the genotypebased scenario. Disease models “Rare”, “Both” and “Common” are described in the section “Population genetics simulation”. Type1 error is set to 5%. Panel 1: 50% of rare variants/haplotypes were assumed to be causal; Panel 2: 20% of rare variants/haplotypes were assumed to be causal; Panel 3: 10% of rare variants/haplotypes were assumed to be causal; Panel 4: empirical type1 error estimate for simulations under the null hypothesis.
Format: PDF Size: 700KB Download file
This file can be viewed with: Adobe Acrobat Reader
Application to central corneal thickness GWAS data set
A total of 552318 SNPs were mapped using the hierarchical mapping algorithm described in the “Methods” section. As a result, we obtained 36146 genes and betweengene blocks. Regions that reached genomewide significance (1.38E6 after Bonferroni correction) for at least one of the four applied tests are presented in Table 2. We identified all of the significant regions reported by Vithana et al. [41]: COL8A2 gene, an interval between the genes RXRA and COL5A1 and a region near ZNF469 gene. As can be seen from Table 2, genotypebased SKAT and MinPval tests achieved genomewide significance for all the five regions listed, whereas SumPval failed to reach genomewide significance for both RXRACOL5A1 and C7orf42 regions. This highlights that MinPval approach performed better than the SumPval approach. Haplotypebased SKAT failed to identify any association signal for four out of the five regions. From our results it is clear that for this data set genotypes were more relevant for identifying associated regions. Judging from the performance of the Haplotype SKAT there is no evidence of association of haplotypes with a phenotype except for COL8A2 gene. It is also of interest to compare our results to those reported by Vithana et al. [41] for singleSNP analysis. For the two out of three regions genotypebased SKAT and MinPval methods outperformed the singleSNP analysis and yielded lower pvalues, which may be explained by the utilization of linkage disequilibrium (LD) within a region to boost power. Given the sensitivity of SumPval approach to underperformance of one of the underlying tests, this method showed higher pvalues. To justify our assumption of bivariate normality for calculation of pvalues for the proposed tests, we used the ShapiroWilk test. The corresponding pvalues for the significant regions are presented in the (Additional file 1: Table S1) in the first row. All the pvalues, except those for COL8A2TRAPPC3 region, are nonsignificant at the genomewide level, which suggests there is no evidence against the bivariate normality assumption. The ShapiroWilk test of COL8A2TRAPPC3 region yielded a marginally significant pvalue on the genomewide level. Hence, the pvalue for this region in the Table 2 is based on permutations.
Table 2. The results of the combined SiMES and SINDI data analysis and the singleSNP pvalues from the original article
Table 3 shows the results of the replication analysis. For the region COL8A2TRAPPC3 the reported replication pvalue is based on permutations, whereas for other regions there was no evidence against bivariate normality assumption (Additional file 1: Table S1, second row). As can be seen from Table 3, only RXRACOL5A1 region was significant after Bonferroni correction for all the tests except for the haplotypebased SKAT. Our replication results are consistent with those of Cornes et al. [54] who found strong evidence of association of multiple SNPs within the RXRACOL5A1 region, and marginal significance of COL8A2 SNP rs96067. It is worth noting that in our analysis C7orf42 gene, which was not identified by Vithana et al. [41], reached genomewide significance in SiMES + SINDI data set and had moderate pvalue in the replication dataset. Cornes et al. [54] found this gene to be significant in a metaanalysis of SiMES, SINDI, 1883 samples from the Singapore Chinese Eye Study and 798 samples from the Beijing Eye Study [55]. The role of C7orf42 gene in central corneal thickness (CCT) phenotype requires further investigation. The results of our analysis suggest that RXRACOL5A1 region may have an impact on CCT phenotype.
Table 3. Replication results on Chinese samples from the Singapore Indian Chinese cohort eye study
In addition to the genebased analysis, we tried to replicate the four genomewide significant SNPs found by Vithana et al. [41] in our Chinese samples using singleSNP analysis. Having tested an association of these SNPs with CCT trait using trend test within a linear additive model adjusting for age, gender and the first ten principal components, we found that none of the SNPs was significant on the corrected type1 error rate 0.0125 = 0.05/4. This result suggests that genebased replication may be a more powerful strategy than singleSNP replication.
In addition to the main genomewide analysis of SiMES + SINDI data set, we applied the proposed methods with a different pair of underlying tests to the three regions reported by Vithana et al. [41]. Both MinPval and SumPval identified the three regions on genomewide significance level. This result suggests that our combined approaches work as well with other underlying tests (for more details, see Additional file 1).
Discussion
When the underlying disease model is unknown, combining statistical tests tailored for different disease scenarios may be a much better strategy than application of a statistical test designed for one specific disease model. In this article we have described the two approaches of combining genotype and haplotypebased statistical tests. The results of theoretical power considerations, population genetics simulations and real data analysis showed strong performance of MinPval approach for different disease scenarios, whereas SumPval method was shown to perform poorly when one of the underlying tests had low power. Our analysis of SiMES + SINDI identified the three regions found by Vithana et al. [41], and additionally, the C7orf42 gene. The replication analysis confirmed an association of RXRACOL5A1 region, which is consistent with the results of Cornes et al. [54], and showed a moderate pvalue for C7orf42 gene. The analysis of real data highlighted the applicability of our combined approaches to real association studies.
In our simulations the Haplotype SKAT was the most powerful test in many cases, but in real data analysis it performed the worst. It is not known beforehand whether a genotype or a haplotypebased test would perform better; hence, our proposal to apply a combined approach is a robust choice. Indeed, MinPval did well in both simulations and real data. This emphasizes the major point of the combined strategy: MinPval may have slightly lower power when a disease model fits Haplotype SKAT and higher power when the disease model is closer to the second underlying tests. One of the possible reasons for the apparent inconsistency of Haplotype SKAT performance may be that for “Rare” and “Both” simulation models we assumed that rare variants bear the major association signal whereas in the real data only common SNPs were present. However, Haplotype SKAT performed well even for “Common” model when a common SNP was causal. We suppose that for this scenario genotype association translated into an association of haplotypes with a phenotype, which is possible if common SNPs within a region are in high LD with each other. On the other hand, if a causal common SNP within a region is in low LD with other common SNPs within a region then under a genotypebased disease scenario haplotypebased test may have much lower power than a genotypebased test which is observed in the results of the real data analysis.
The methods proposed in this study may be easily generalized to multiple statistical tests, namely, instead of two underlying tests it is possible to apply more tests and combine all of them via the described methodology. In this case the arguments for theoretical pvalue calculation for the proposed approaches can be extended in a straightforward manner.
Recently Derkach et al. [56] investigated the performance of the combined approaches, namely, the minimum of pvalues and the Fisher pvalue combination, for rare variants association scenarios. Although the approaches we propose are similar, our major idea is different. We combine two test statistics for the purpose of widening the set of alternatives for which our test is powerful; thus, we choose the underlying tests designed for very different phenotype models, whereas Derkach et al. [56] used linear and quadratic tests which are likely to be both powerful under many models. As a result, our conclusions are different from those of Derkach et al. [56]. For example, the authors stated that “hybrid test statistics provide much needed robustness in terms of power for association tests”, whereas we observed that only minimum pvalue approach really preserves power when one of the underlying tests underperforms. Secondly, the authors found that in many cases Fisher method outperforms both of the underlying tests, and the minimum pvalue approach. However, from our work it is clear that SumPval (which is similar to the Fisher pvalue combination) outperforms all the three tests only when both of the underlying tests have comparable power which is unlikely if the two underlying tests are deliberately chosen to fit very different phenotype models.
One of the limitations of the proposed approaches is the need to use permutations. For theoretical pvalue calculation both SumPval and MinPval require a correlation coefficient to be estimated via permutations. Moreover, permutations need to be applied when asymptotic distributions of the underlying test statistics are unknown or inadequate to describe the empirical distributions.
The described methodologies may be extended to preserve power under other disease models. For example, the combination of rarevariants and commonvariants statistical tests applied to a sequenced region may preserve high power when either only rare or only common variants are associated with a phenotype. However, it is not known how the combined approaches will perform if both common and rare variants are associated with phenotype.
Conclusions
In this study we have investigated the performance of combined haplotype and genotypebased tests for the purpose of preserving high power under both genotype and haplotype disease scenarios. Based on theoretical power calculations, population genetics simulations and analysis of the real data set we have illustrated high performance and potential utility of combined approaches for association studies.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
SZ, AS and AT conceived the study. SZ and AT designed the experiments. SZ conducted the experiments and performed the analysis. TYW, TA, ENV, and CCK provided the GWAS data. SZ and AT wrote the manuscript. SZ, TYW, TA, EV, KCC, AS and AT approved the manuscript.
Funding
This work was supported by the Agency for Science, Technology and Research (A*STAR), Singapore. The first author is a recipient of the Singapore International Graduate Award.
References

Green EK, Grozeva D, Sims R, Raybould R, Forty L, GordonSmith K, Russell E, St. Clair D, Young AH, Ferrier IN, et al.: DISC1 exon 11 rare variants found more commonly in schizoaffective spectrum cases than controls.
Am J Med Genet B Neuropsychiatr Genet 2011, 156(4):490492. Publisher Full Text

Norton N, Li D, Rieder Mark J, Siegfried Jill D, Rampersaud E, Züchner S, Mangos S, GonzalezQuintana J, Wang L, McGee S, et al.: Genomewide studies of copy number variation and exome sequencing identify rare variants in BAG3 as a cause of dilated cardiomyopathy.
The American Journal of Human Genetics 2011, 88(3):273282. Publisher Full Text

Ramagopalan SV, Dyment DA, Cader MZ, Morrison KM, Disanto G, Morahan JM, BerlangaTaylor AJ, Handel A, De Luca GC, Sadovnick AD, et al.: Rare variants in the CYP27B1 gene are associated with multiple sclerosis.
Ann Neurol 2011, 70(6):881886. PubMed Abstract  Publisher Full Text

Xie P, Kranzler HR, Krauthammer M, Cosgrove KP, Oslin D, Anton RF, Farrer LA, Picciotto MR, Krystal JH, Zhao H, et al.: Rare nonsynonymous variants in alpha4 Nicotinic Acetylcholine receptor gene protect against nicotine dependence.
Biol Psychiatry 2011, 70(6):528536. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wu MC, Kraft P, Epstein MP, Taylor DM, Chanock SJ, Hunter DJ, Lin X: Powerful SNPSet analysis for case–control genomewide association studies.
Am J Hum Genet 2010, 86(6):929942. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Caporaso N, Gu F, Chatterjee N, ShengChih J, Yu K, Yeager M, Chen C, Jacobs K, Wheeler W, Landi MT, et al.: Genomewide and candidate gene association study of cigarette smoking behaviors.
PLoS ONE 2009, 4(2):e4653. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Hong MG, Reynolds CA, Feldman AL, Kallin M, Lambert JC, Amouyel P, Ingelsson E, Pedersen NL, Prince JA: Genomewide and genebased association implicates FRMD6 in alzheimer disease.
Hum Mutat 2012, 33(3):521529. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Akey J, Jin L, Xiong M: Haplotypes vs single marker linkage disequilibrium tests: what do we gain?
Eur J Hum Genet 2001, 9:291300. PubMed Abstract  Publisher Full Text

Schaid DJ: Power and sample size for testing associations of haplotypes with complex traits.
Ann Hum Genet 2006, 70(1):116130. PubMed Abstract  Publisher Full Text

Browning SR, Browning BL: Rapid and accurate Haplotype phasing and missingdata inference for wholegenome association studies by use of localized Haplotype clustering.
The American Journal of Human Genetics 2007, 81(5):10841097. Publisher Full Text

Howie BN, Donnelly P, Marchini J: A flexible and accurate genotype imputation method for the next generation of genomewide association studies.
PLoS Genet 2009, 5(6):e1000529. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stephens M, Smith NJ, Donnelly P: A New statistical method for haplotype reconstruction from population data.
Am J Hum Genet 2001, 68(4):978989. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bansal V, Halpern AL, Axelrod N, Bafna V: An MCMC algorithm for haplotype assembly from wholegenome sequence data.
Genome Res 2008, 18(8):13361346. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Bansal V, Libiger O, Torkamani A, Shork JN: Statistical analysis strategies for association studies involving rare variants.

Schatz MC, Delcher AL, Salzberg SL: Assembly of large genomes using secondgeneration sequencing.
Genome Res 2010, 20(9):11651173. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Begnini A, Tessari G, Turco A, Malerba G, Naldi L, Gotti E, Boschiero L, Forni A, Rugiu C, Piaserico S, et al.: PTCH1 gene haplotype association with basal cell carcinoma after transplantation.
Brit J Dermatol 2010, 163(2):364370. Publisher Full Text

DIEUDE P, DAWIDOWICZ K, GUEDJ M, LEGRAIN Y, WIPFF J, HACHULLA E, DIOT E, SIBILIA J, MOUTHON L, CABANE J, et al.: PhenotypeHaplotype Correlation of IRF5 in Systemic Sclerosis: role of 2 Haplotypes in Disease Severity.
J Rheumatol 2010, 37(5):987992. PubMed Abstract  Publisher Full Text

Lambert JC, GrenierBoley B, Harold D, Zelenika D, Chouraki V, Kamatani Y, Sleegers K, Ikram MA, Hiltunen M, Reitz C, et al.: Genomewide haplotype association study identifies the FRMD4A gene as a risk locus for Alzheimer's disease.
Mol Psychiatry 2013, 18(4):461470. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Tregouet DA, Konig IR, Erdmann J, Munteanu A, Braund PS, Hall AS, Groszhennig A, LinselNitschke P, Perret C, DeSuremain M, et al.: Genomewide haplotype association study identifies the SLC22A3LPAL2LPA gene cluster as a risk locus for coronary artery disease.
Nat Genet 2009, 41(3):283285. PubMed Abstract  Publisher Full Text

Bansal V, Bafna V: HapCUT: an efficient and accurate algorithm for the haplotype assembly problem.
Bioinformatics 2008, 24(16):i153i159. PubMed Abstract  Publisher Full Text

Gauderman WJ, Murcray C, Gilliland F, Conti DV: Testing association between disease and multiple SNPs in a candidate gene.
Genet Epidemiol 2007, 31(5):383395. PubMed Abstract  Publisher Full Text

Li M, Fu W, Lu Q: An aggregating UTest for a genetic association study of quantitative traits.
BMC Proceedings 2011, 5(Suppl 9):S23. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Li YM, Xiang Y, Sun ZQ: An entropybased measure for QTL mapping using extreme samples of population.
Hum Hered 2008, 65(3):121128. PubMed Abstract  Publisher Full Text

Bansal V, Libiger O, Torkamani A, Schork N: Statistical analysis strategies for association studies involving rare variants.

IonitaLaza I, Buxbaum JD, Laird NM, Lange C: A New testing strategy to identify rare variants with either risk or protective effect on disease.
PLoS Genet 2011, 7(2):e1001289. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Madsen BE, Browning SR: A groupwise association test for rare mutations using a weighted Sum statistic.
PLoS Genet 2009, 5(2):e1000384. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Neale BM, Rivas MA, Voight BF, Altshuler D, Devlin B, OrhoMelander M, Kathiresan S, Purcell SM, Roeder K, Daly MJ: Testing for an unusual distribution of rare variants.
PLoS Genet 2011, 7(3):e1001322. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Price AL, Kryukov GV, de Bakker PIW, Purcell SM, Staples J, Wei LJ, Sunyaev SR: Pooled association tests for rare variants in exonresequencing studies.
Am J Hum Genet 2010, 86(6):832838. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Jin L, Zhu W, Guo J: Genomewide association studies using haplotype clustering with a new haplotype similarity.
Genet Epidemiol 2010, 34(6):633641. PubMed Abstract  Publisher Full Text

Sha Q, Dong J, Jiang R, Zhang S: Tests of association between quantitative traits and haplotypes in a reduceddimensional space.
Ann Hum Genet 2005, 69(6):715732. PubMed Abstract  Publisher Full Text

Zaykin DV, Westfall PH, Young SS, Karnoub MA, Wagner MJ, Ehm MG: Testing association of statistically inferred haplotypes with discrete and continuous traits in samples of unrelated individuals.
Hum Hered 2002, 53(2):7991. PubMed Abstract  Publisher Full Text

Guo W, Lin S: Generalized linear modeling with regularization for detecting common disease rare haplotype association.
Genet Epidemiol 2009, 33(4):308316. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Li Y, Byrnes AE, Li M: To identify associations with rare variants, just WHaIT: weighted haplotype and imputationbased tests.
Am J Hum Genet 2010, 87(5):728735. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Zhu X, Feng T, Li Y, Lu Q, Elston RC: Detecting rare variants for complex traits using family and unrelated data.
Genet Epidemiol 2010, 34(2):171187. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Stouffer SA, Suchman EA, Devinney LC, Star SA, Williams RMJ: Adjustment during army life. In The American soldier. 1st edition. Princeton, NJ: Princeton Univ; 1949.

Tippet LHC: The method of statistics. London: Williams and Northgate; 1931.

King CR, Rathouz PJ, Nicolae DL: An evolutionary framework for association testing in resequencing studies.
PLoS Genet 2010, 6(11):e1001202. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Boyko AR, Williamson SH, Indap AR, Degenhardt JD, Hernandez RD, Lohmueller KE, Adams MD, Schmidt S, Sninsky JJ, Sunyaev SR, et al.: Assessing the evolutionary impact of amino acid mutations in the human genome.
PLoS Genet 2008, 4(5):e1000083. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Gutenkunst RN, Hernandez RD, Williamson SH, Bustamante CD: Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data.
PLoS Genet 2009, 5(10):e1000695. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wu Michael C, Lee S, Cai T, Li Y, Boehnke M, Lin X: Rarevariant association testing for sequencing data with the sequence kernel association test.
Am J Hum Genet 2011, 89(1):8293. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Vithana EN, Aung T, Khor CC, Cornes BK, Tay WT, Sim X, Lavanya R, Wu R, Zheng Y, Hibberd ML, et al.: Collagenrelated genes influence the glaucoma risk factor, central corneal thickness.
Hum Mol Genet 2011, 20(4):649658. PubMed Abstract  Publisher Full Text

Lavanya R, Jeganathan VSE, Zheng Y, Raju P, Cheung N, Tai ES, Wang JJ, Lamoureux E, Mitchell P, Young TL, et al.: Methodology of the Singapore Indian Chinese cohort (SICC) Eye study: quantifying ethnic variations in the epidemiology of eye diseases in Asians.
Ophthalmic Epidemiol 2009, 16(6):325336. PubMed Abstract  Publisher Full Text

Foong AWP, Saw SM, Loo JL, Shen S, Loon SC, Rosman M, Aung T, Tan DTH, Tai ES, Wong TY: Rationale and methodology for a populationbased study of Eye diseases in Malay people: the Singapore Malay Eye study (SiMES).
Ophthalmic Epidemiol 2007, 14(1):2535. PubMed Abstract  Publisher Full Text

Su DHW, Wong TY, Foster PJ, Tay WT, Saw SM, Aung T: Central corneal thickness and its associations with ocular and systemic factors: the Singapore Malay Eye study.
Am J Ophthalmol 2009, 147(4):709716.e701. PubMed Abstract  Publisher Full Text

Wong TCE, et al.: Prevalence and causes of low vision and blindness in an urban malay population: the singapore malay eye study.
Arch Ophthalmol 2008, 126(8):10911099. PubMed Abstract  Publisher Full Text

Zhao J, Gupta S, Seielstad M, Liu J, Thalamuthu A: Pathwaybased analysis using reduced gene subsets in genomewide association studies.
BMC Bioinformatics 2011, 12:17. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D: Principal components analysis corrects for stratification in genomewide association studies.
Nat Genet 2006, 38(8):904909. PubMed Abstract  Publisher Full Text

Thalamuthu A, Zhao J, Keong G, Kondragunta V, Mukhopadhyay I: Association tests for rare and common variants based on genotypic and phenotypic measures of similarity between individuals.
BMC Proceedings 2011, 5(Suppl 9):S89. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Freedman D, Lane D: A nonstochastic interpretation of reported significance levels.
Journal of Business and Economic Statistics 1983, 1:292298.

Ter Braak CJF: Permutation versus bootstrap significance tests in multiple regression and ANOVA. In Bootstrapping and related techniques. Edited by Jockel KH, Rothe G, Sendler W. Berlin: Springer; 1992.

Anderson MJ, Legendre P: An empirical comparison of permutation methods for tests of partial regression coefficients in a linear model.
J Stat Comput Sim 1999, 62(3):271303. Publisher Full Text

Wagner BD, Zerbe GO, Mexal S, Leonard SS: Permutationbased adjustments for the significance of partial regression coefficients in microarray data analysis.
Genet Epidemiol 2008, 32(1):18. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Cornes BK, Khor CC, Nongpiur ME, Xu L, Tay WT, Zheng Y, Lavanya R, Li Y, Wu R, Sim X, et al.: Identification of four novel variants that influence central corneal thickness in multiethnic Asian populations.
Hum Mol Genet 2012, 21(2):437445. PubMed Abstract  Publisher Full Text

Zhang H, Xu L, Chen C, Jonas J: Central corneal thickness in adult Chinese. Association with ocular and general parameters. The Beijing Eye Study.
Graefe’s Archive for Clinical and Experimental Ophthalmology 2008, 246(4):587592. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Derkach A, Lawless JF, Sun L: Robust and powerful tests for rare variants using Fisher’s method to combine evidence of association from two or more complementary tests.
Genet Epidemiol 2013, 37(1):110121. PubMed Abstract  Publisher Full Text