Email updates

Keep up to date with the latest news and content from BMC Genetics and BioMed Central.

Open Access Research article

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

Ludwig A Hothorn12*, Ondrej Libiger2 and Daniel Gerhard1

Author affiliations

1 Institute of Biostatistics, Leibniz University Hannover, D–30419 Hannover, Germany

2 Scripps Genomic Medicine, The Scripps Research Institute, La Jolla, CA–92037, USA

For all author emails, please log on.

Citation and License

BMC Genetics 2012, 13:59  doi:10.1186/1471-2156-13-59

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2156/13/59


Received:4 August 2011
Accepted:18 July 2012
Published:18 July 2012

© 2012 Hothorn et al.; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

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.

Keywords:
Genetic association study; Quantitative traits; Interaction; Variance heterogeneity

Author’s summary

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.

Background

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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M1">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M2">View MathML</a> against the alternative that a difference exists between any pair of variances: <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M3">View MathML</a>. The test statistic consists of a quadratic form

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M4">View MathML</a>

(1)

(with J=3) using the robust Levene residuals

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M5">View MathML</a>

(2)

[8] with nj quantitative trait observations Yij per genotype j. The <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M6">View MathML</a> is F-distributed with df1=Jand df2=NJ.

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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M7">View MathML</a> vs. <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M8">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M9">View MathML</a> vs. <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M10">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M11">View MathML</a> vs. <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M12">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M13">View MathML</a> denotes the variance estimator for the pooled groups jj. 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].

Methods

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:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M14">View MathML</a>

(3)

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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M15">View MathML</a> vs. <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M16">View MathML</a><a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M17">View MathML</a> vs. <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M18">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M19">View MathML</a> vs. <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M20">View MathML</a>[10], are formulated in the 3×3 contrast matrix:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M21">View MathML</a>

(4)

A priori it is unknown which elementary test is mostly under the alternative. Therefore, the maximum of the test statistics <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M22">View MathML</a> 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 ν=<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M23">View MathML</a> degrees of freedom and a correlation matrix R=(ρkk) given by its elements

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M24">View MathML</a>

(5)

Under the above assumptions the correlations ρkkdepend 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M25">View MathML</a>[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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M26">View MathML</a> 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 MCTPairs can be used by means of the following 3×3 contrast matrix:

<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M27">View MathML</a>

(6)

Results and discussion

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 MCTPairsand MCTAveapproaches allows elementary decisions additionally. The related simulations for the additive and recessive modes are presented in the Appendix.

Table 1. Size and power comparison of the Levene test and two contrast alternatives given an dominant mode of inheritance

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 <10−8were 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.

thumbnailFigure 1. Boxplots including raw data for two selected SNPs and phenotypes of the Bogalusa Heart Study.

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 MCTAve) 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 (<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M28">View MathML</a>), and the pooled variance for the two remaining genotype groups associated with the presence of at least one non-risk allele (<a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M29">View MathML</a>). 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 <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M30">View MathML</a> was significantly higher compared to homozygotes <a onClick="popup('http://www.biomedcentral.com/1471-2156/13/59/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1471-2156/13/59/mathml/M31">View MathML</a>. 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 MCTPairs) 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.

Table 2. P values for original and multiple contrast Levene-type tests

thumbnailFigure 2. QQ- and Manhattan-plot for the phenotype waist circumference in the Bogalusa Heart Study.

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].

Additional file 1. The R code.

Format: PDF Size: 8KB Download file

This file can be viewed with: Adobe Acrobat ReaderOpen Data

Additional file 2. The data settg.rdacontains the subject-specific diastolic blood pressure data and the related genotype levels of the SNP rs12607553.

Format: RDA Size: 1KB Download fileOpen Data

Conclusions

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.

Appendix

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

Table 3. Size and power comparison of the Levene test and two contrast alternatives given an additive mode of inheritance

Table 4. Size and power comparison of the Levene test and two contrast alternatives given an recessive mode of inheritance

Authors’ contributions

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.

Acknowledgements

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.

References

  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. OpenURL

  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. OpenURL

  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. OpenURL

  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. OpenURL

  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. OpenURL

  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. OpenURL

  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. OpenURL

  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. OpenURL

  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. OpenURL

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

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

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

    J Qual Technol 2009, 41:60-65. OpenURL

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

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

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

    Ann Stat 1990, 18:465-468. OpenURL

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

    R News 2002, 2:14-17. OpenURL

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