Model-specific tests on variance heterogeneity for detection of potentially interacting genetic loci

Background

Trait variances among genotype groups at a locus are expected to differ in the presence of an interaction between this locus and another locus or environment. A simple maximum test on variance heterogeneity can thus be used to identify potentially interacting single nucleotide polymorphisms (SNPs).

Results

We propose a multiple contrast test for variance heterogeneity that compares the mean of Levene residuals for each genotype group with their average as an alternative to a global Levene test. We applied this test to a Bogalusa Heart Study dataset to screen for potentially interacting SNPs across the whole genome that influence a number of quantitative traits. A user-friendly implementation of this method is available in the R statistical software package multcomp.

Conclusions

We show that the proposed multiple contrast test of model-specific variance heterogeneity can be used to test for potential interactions between SNPs and unknown alleles, loci or covariates and provide valuable additional information compared with traditional tests. Although the test is statistically valid for severely unbalanced designs, care is needed in interpreting the results at loci with low allele frequencies.

Interactions among alleles at variant sites in the genome or between alleles and the environment likely play an important role in determining complex traits such as blood pressure. However, sets of interacting loci are difficult to identify due to the large number of potential interactions that need to be tested. One approach that circumvents this difficulty is to identify loci that appear to take part in an interaction although their partners with which they interact are unknown. A SNP locus containing an allele that interacts with other alleles or the environment can be identified by the existence of a statistically significant difference in the variance of quantitative trait values among individuals who possess zero, one or two alleles at the locus. We describe an extension of Levene’s test, which was proposed to test variance heterogeneity. This new test has the advantage of providing information regarding the effect of specific alleles on variance heterogeneity, which can lead to formulating concrete, biologically relevant hypotheses about interacting alleles rather than just loci while controlling for type I error rate.

Statistical association between a biallelic marker and a quantitative trait is usually tested using either a two degree of freedom F-test in the one-way analysis of variance (ANOVA) [1], or a one degree of freedom F-test in a linear regression [2]. These approaches compare the means of quantitative trait values at genotype categories associated with a SNP locus (i.e., homozygous for major allele, heterozygous and homozygous for minor allele). While ANOVA is sensitive to any global heterogeneity, linear regression test is sensitive to the presence of an additive mode of inheritance. Less attention has been given to comparing the variances in the quantitative trait values associated with different genotype categories. Recently, [3] proposed using a standard Levene test [4] to identify variance heterogeneity due to potential interaction between a given locus and another allele at the same locus, alleles at different loci or the environment. They compared three global tests, namely the Bartlett-test, a rank modification of Bartlett test and Levene test particularly for non-normal distributed variables. Differences among the variances of quantitative trait values at each genotype category (denoted with j=0,1,2 interacting alleles) may reflect an interaction [5]. In contrast to approaches that explicitly test specific gene-gene, e.g. by Bayesian partition methods [6] or gene-environment interactions, e.g. by multiple regression methods [7], methods that assess variance heterogeneity can be used to uncover loci that are not previously known to interact.

Levene test [8] tests a global null hypothesis against the alternative that a difference exists between any pair of variances: . The test statistic consists of a quadratic form

(with J=3) using the robust Levene residuals

[8] with _nj quantitative trait observations _Yij per genotype j. The is F-distributed with d_f1=Jand d_f2=N− J.

This test is known to be relatively robust when data are not normally distributed. However, the main disadvantage of Levene’s test is that it can only be used to determine whether the group-specific variances differ among each other. In order to obtain a biologically or clinically relevant interpretation of the results, it is often valuable to additionally determine which pairs of genotype categories in particular exhibit statistically significant variance heterogeneity.

To this end, [3] considered using three two-sample df−1 tests for the three comparisons vs. vs. , and vs. , where denotes the variance estimator for the pooled groups j^j″. However, these multiple tests do not control for the family-wise type I error rate α.

In this paper, we propose a Levene-type multiple contrast test, a novel approach comprised of a global test on variance heterogeneity as well as the three specific tests on pairwise variance heterogeneity using a maximum test of linear forms. We apply this test in a genome-wide fashion using a Bogalusa Heart Study dataset [9].

A Levene-type multiple contrast test

For the Levene-type transformed variable _Zijand the factor genotype with the levels j=0,1,2 the following contrast test for the one-way layout is used:

where S is the root of the common mean square error estimate of the linear model on the transformed data _Zji (absolute Levene residuals). Contrast coefficients _ckjare used for the k=3 comparisons between the means of Levene residuals for each genotype with the average of the means. These comparisons coding for vs. vs. , and vs. [10], are formulated in the 3×3 contrast matrix:

A priori it is unknown which elementary test is mostly under the alternative. Therefore, the maximum of the test statistics is used, implying a family-wise type-I-error rate αfor all of the three comparisons.

Under the approximate assumption of multivariate normal distributed errors with a homogeneous global variance of the transformed variable _Zij, the vector ^(_Tj)″ follows jointly a tri-variate t-distribution with ν= degrees of freedom and a correlation matrix R=(_ρk^k″) given by its elements

Under the above assumptions the correlations _ρk^k″depend only on the contrast coefficients _ckjand the sample sizes _nj. This approach controls the familywise error rate αand reveals a reasonable power for unbalanced designs as long as the above assumptions hold true. It provides a global decision whenever any of the contrasts is under the alternative and additionally the elementary decisions by multiplicity-adjusted p-values for the three specific comparisons. Although simultaneous confidence intervals for both differences and ratios to the average are available as well [11], they will not be recommended since a genetic interpretation for the transformed variable _Zji is difficult. Recently, simultaneous confidence intervals for the pairwise ratios of variances were proposed using a maximum test on jackknifed [12]. This approach would provide an alternative when modified for arbitrarily unbalanced designs.

The question arises whether the quadratic form of the Levene test or the maximum contrast version of linear forms is more powerful. A general answer will not be available. For the comparison of normally distributed means, the least favorable configuration approach for reveals higher power for the maximum of linear contrasts relative to the quadratic form of the F-test of the one-way analysis of variance [13]. Alternatively, a multiple contrast test based on pairwise comparisons of the genotype-specific variances MC^TPairs can be used by means of the following 3×3 contrast matrix:

Simulation Study

To validate the type-I-error control and compare the power between the multiple contrast test on the Levene scores and the Levene test, a small simulation study was conducted. Observations are sampled from a standard normal distribution for the three genotype groups. Overall sample size N is varied from 25 to 100, allocating the observations under assumption of Hardy-Weinberg equilibrium with allele frequencies p for the recessive allele of 0.5 and 0.75. To generate differences in variances, the variance is increased with the number of dominant alleles by a value of δ for an additive mode of inheritance (_σ0=1,_σ1=1 + δ,_σ2=1 + 2δ). For a recessive and dominant model, the value of δcorresponding to the intermediate genotype is changed to 0 or 2 δ, respectively. For each parameter setting 10,000 simulation replications are performed. In Table 1 for the three approaches both control of type-I-error and comparable powers from a global test perspective can be concluded for the dominant mode of inheritance. However, the new MC^TPairsand MC^TAveapproaches allows elementary decisions additionally. The related simulations for the additive and recessive modes are presented in the Appendix.

Evaluation of a real data example

We evaluated the utility of the proposed test using data from Bogalusa Heart Study (BHS). The longitudinal study included genotype information on 525 unrelated individuals of European descent at 545,821 SNPs [12]. Twelve clinically-relevant quantitative traits were also measured for each study participant. We applied the multiple contrast test to each SNP throughout the genome. A number of SNP loci with a p-value smaller than <1⁰⁻⁸were identified for several quantitative traits whereas the number of SNPs with p-values below this threshold varied for different quantitative traits. For example, 23 SNP loci yielded p-values lower than that threshold when considering heart rate, while only 5 SNPs reached similar significance in the analysis involving body mass index. A thorough examination of the results showed that several significant tests occurred when the risk allele homozygote genotype group contained only a handful, i.e. 1,2,3, data points. Although the approach described in this paper is statistically valid in these situations, we selected only tests in which each genotype group contained at least one percent of all data points for further study. In the following text, we showcase two instances of variance heterogeneity in greater detail: rs3760124 and waist circumference, as well as rs12607553 and diastolic blood pressure; see the related box-plots in Figure 1.

The boxplots show relatively symmetric distributions of quantitative trait values in all genotype groups, thus ruling out the presence of outliers or extremely skewed distributions of trait values as sources of the observed variance heterogeneity. Furthermore, all three genotype categories contain a relatively large number of observed trait values resulting in reliable variance estimates.

It can be seen from Table 2 that our proposed test (labeled as MC^TAve) yielded a result with greater statistical significance compared with the Levene global test. More importantly however, it allowed us to interpret the effect of the SNP alleles on variance heterogeneity. When we considered waist circumference, we found the variance heterogeneity to be most significant between the risk allele homozygous genotype group (), and the pooled variance for the two remaining genotype groups associated with the presence of at least one non-risk allele (). In addition, we can see that the variance associated with the risk allele homozygotes is significantly higher in this example. Conversely, in the second example involving diastolic blood pressure, the variance in the heterozygotes was significantly higher compared to homozygotes . Such insights may prove to be valuable in forming hypotheses for further research, and can be obtained with the Levene-type multiple contrast test that we propose here. Furthermore, pairwise comparisons (labeled as MC^TPairs) can also be performed (see Table 2). From the genome-wise analysis of the endpoint waist circumference the quantile-quantile plot and the Manhattan plot are presented in Figure 2 whereas SNPs with 5 or more individuals who are homozygous for the minor allele were considered. We see the expected inflationary behavior on a genome-wise level and the relatedness of the SNPs within the blocks. But only one SNP reveals is significant interaction on a global Bonferroni level.

An example R code for testing variance heterogeneity at a single SNP is provided in the Additional files 1 and 2. The multiplicity-adjusted p-values can be estimated by means of the R package multcomp [14]. Alternatively, a SAS procedure GLIMMIX can be used for a resampling based estimation of multiplicity-adjusted p-values [15].

The important issue of missing heritability, which refers to the fact that common SNPs identified by genome-wide association studies as associated with a disease collectively explain only a small portion of the prevalence of this disease, may be due, in part, to the presence of unknown interactions among alleles at various SNP loci or environment, that affect the disease. The identification of such interactions is difficult, primarily because of the large number of potentially interacting pairs, trios, etc. of alleles and environmental variables that need to be tested. A feasible alternative, as suggested by [5], is to test individual loci for the evidence of their involvement in an interaction with other alleles, loci or covariates. The idea of assessing variance heterogeneity between three genotype groups at a particular SNP locus as evidence of a potential interaction is appealing for its simplicity. The Levene-type maximum contrast test proposed in this paper allows one to not only test for global variance heterogeneity, but also perform groupwise test that allows one to elucidate the effect of the individual alleles on quantitative trait variance. While this is an advantage over the standard Levene test, the price to pay is increased computing time. However, we were able to perform a genome-wide analysis of variance heterogeneity involving >500 individuals in a matter of minutes on a laptop computer. Parallelization can also be used to substantially decrease computation time requirements. R code implementing this test is available as part of the multcomp package.

Even the analysis of real data example illustrates the low specificity of the identified potential interactions. Care needs to be exercised in interpreting the results of this test in cases of low frequency variants or missing trait data, when one or more of the genotype groups contains an extremely small number of observed trait values. The issues surrounding the sensitivity and specificity of this approach in these potentially common cases is an area that needs further work.

The simulation results for the additive and recessive mode of inheritance are presented in Tables 3 and 4.

LAH and DG derived the method. LAH and OL developed the software. OL performed the GWA analysis. LAH, OL and DG wrote the manuscript. All authors read and approved the final manuscript.

We appreciate the two anonymous reviewers for their insightful and constructive comment. The work for the first author was partly supported by the German Science Foundation grand DfG-HO1687.

1. Liu YZ, Pei YF, Guo YF, Wang L, Liu XG, Yan H, Xiong DH, Zhang YP, Levy S, Li J, Haddock CK, Papasian CJ, Xu Q, Ma JZ, Payne TJ, Recker RR, Li MD, Deng HW: Genome-wide association analyses suggested a novel mechanism for smoking behavior regulated by IL15.

   Mol Psychiatry 2009, 14(7):668-680.

2. Zhao JH, Li MY, Bradfield JP, Zhang HT, Mentch FD, Wang K, Sleiman PM, Kim CE, Glessner JT, Hou CP, Keating BJ, Thomas KA, Garris ML, Deliard S, Frackelton EC, Otieno FG, Chiavacci RM, Berkowitz RI, Hakonarson H, Grant SFA: The role of height-associated loci identified in genome wide association studies in the determination of pediatric stature.

   Bmc Med Genet 2010, 11:96.

3. Struchalin MV, Dehghan A, Witteman JCM, van Duijn, Aulchenko YS: Variance heterogeneity analysis for detection of potentially interacting genetic loci: method and its limitations.

   Bmc Genet 2010, 11:92.

4. Levene H: Robust tests for equality of variances. In Contributions to Probability and Statistics. Edited by Olkin I. Palo Alto, CA: Stanford University Press; 1960:278-292.

5. Pare G, Cook NR, Ridker PM, Chasman DI: On the use of variance per genotype as a tool to identify quantitative trait interaction effects: A report from the women’s genome health study.

   Plos Genet 2010, 6(6):e1000981.

6. Zhang Y, Jiang B, Zhu J, Liu JS: Bayesian models for detecting epistatic interactions from genetic data.

   Ann Human Genet 2011, 75(Part 1):183-193.

7. Erhardt V, Bogdan M, Czado C: Locating multiple interacting quantitative trait loci with the zero-inflated generalized poisson regression.

   Stat Appl Genet Mol Biol 2010, 9:26.

8. Gastwirth JL, Gel YR, Miao WW: The impact of levene’s test of equality of variances on statistical theory and practice.

   Stat Sci 2009, 24(3):343-360.

9. Smith EN, Chen W, Kahonen M, et al.: Longitudinal Genome-Wide Association of Cardiovascular Disease Risk Factors in the Bogalusa Heart Study.

   Plos Genet 2010, 6(9):e1001094.

10. Hasler M, Hothorn LA: Multiple Contrast Tests in the Presence of Heteroscedasticity.

    Biometrical J 2008, 50(5):793-800.

11. Djira GD, Hothorn LA: Detecting Relative Changes in Multiple Comparisons with an Overall Mean.

    J Qual Technol 2009, 41:60-65.

12. Rublik F: A note on simultaneous confidence intervals for ratio of variances.

    Commun Stat-Theory Methods 2010, 39(6):1038-1045.

13. Hayter AJ, Liu W: The Power Function of the Studentised Range Test.

    Ann Stat 1990, 18:465-468.

14. Bretz F, Hothorn T, Westfall P: On multiple comparisons in R.

    R News 2002, 2:14-17.

15. SAS Institute Inc. 2008: SAS/STAT® 9.2, User′s Guide. Cary NC: SAS Institute Inc.;