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

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.

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 ${\sigma }_j^2$ 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 $H_0\text{:}{\sigma }_0^2={\sigma }_1^2={\sigma }_2^2$ against the alternative that a difference exists between any pair of variances: $H_1\text{:}\exists j,j^{″}\text{:}{\sigma }_j^2\ne {\sigma }_{j^{″}}^2,j\ne j^{″}$. The test statistic consists of a quadratic form

$$T_{\mathrm{Levene}}^2=\frac{(N-J)\sum _{j=0}^J{n}_j{(Z_{\mathrm{j.}}-\mathrm{Z..})}^2}{(J-1)\sum _{j=0}^J{(Z_{\mathrm{ji}}-Z_{\mathrm{j.}})}^2},N=\sum _{j=0}^J{n}_j$$

(1)

(with J=3) using the robust Levene residuals

$$Z_{\mathrm{ji}}=\mathrm{abs}(Y_{\mathrm{ji}}-\mathrm{Median}(Y_{\mathrm{j.}}\left)\right),$$

(2)

[8] with _{n j} quantitative trait observations _{Y ij} per genotype j. The $T_{\mathrm{Levene}}^2$ 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 ${\sigma }_0^2$ vs. ${\sigma }_{12}^2$${\sigma }_1^2$ vs. ${\sigma }_{02}^2$, and ${\sigma }_2^2$ vs. ${\sigma }_{01}^2$, where ${\sigma }_{j{j}^{″}}^2$ 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].

Methods

A Levene-type multiple contrast test

For the Levene-type transformed variable _{Z ij} and the factor genotype with the levels j=0,1,2 the following contrast test for the one-way layout is used:

$$T_k^{\mathrm{Levene}}=\frac{\sum _{j=1}^J{c}_{\mathrm{kj}}{\bar{Z}}_{\mathrm{j.}}}{S\sqrt{\sum _{j=1}^J{c}_{\mathrm{kj}}^2/n_j}}$$

(3)

where S is the root of the common mean square error estimate of the linear model on the transformed data _{Z ji} (absolute Levene residuals). Contrast coefficients _{c kj} are 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 ${\sigma }_0^2$ vs. ${\sigma }_{12}^2$${\sigma }_1^2$ vs. ${\sigma }_{02}^2$, and ${\sigma }_2^2$ vs. ${\sigma }_{01}^2$[10], are formulated in the 3×3 contrast matrix:

$$\left(c_{\mathrm{kj}}\right)=\left(\begin{array}{lll}-1& 0.5& 0.5\\ 0.5& -1& 0.5\\ 0.5& 0.5& -1\end{array}\right).$$

(4)

A priori it is unknown which elementary test is mostly under the alternative. Therefore, the maximum of the test statistics $max\left(T_j\right)$ 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 Z _ij , the vector (T _j )^″ follows jointly a tri-variate t-distribution with ν=$\sum _{j=1}^J(n_j-1)$ degrees of freedom and a correlation matrix R=( _{ρ k k} ^″) given by its elements

$${\rho }_{k{k}^{″}}=\frac{\sum _{j=1}^J{c}_{\mathrm{kj}}{c}_{k^{″}j}/n_j}{\sqrt{\left(\sum _{j=1}^J{c}_{\mathrm{kj}}^2/n_j\right)\left(\sum _{j=1}^J{c}_{k^{″}j}^2/n_j\right)}}.$$

(5)

Under the above assumptions the correlations _{ρ k k} ^″depend only on the contrast coefficients _{c kj} and the sample sizes _{n j} . 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 _{Z ji} is difficult. Recently, simultaneous confidence intervals for the pairwise ratios of variances were proposed using a maximum test on jackknifed $\log \left(s_j^2\right)$[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 ${\mu }_j-\overline{\mathrm{\mu .}}$ 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:

$$\left(c_{\mathrm{kj}}\right)=\left(\begin{array}{lll}-1& 1& 0\\ 0& -1& 1\\ -1& 0& 1\end{array}\right).$$

(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 MC^TPairsand MC^TAveapproaches 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 <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.

Figure 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 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 (${\sigma }_{\mathrm{TT}}^2$), and the pooled variance for the two remaining genotype groups associated with the presence of at least one non-risk allele (${\sigma }_{\mathrm{CC},\mathrm{CT}}^2$). 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 ${\sigma }_{\mathrm{AG}}^2$ was significantly higher compared to homozygotes ${\sigma }_{\mathrm{AA},\mathrm{GG}}^2$. 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.

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

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

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