Abstract
Introduction: The standard way of using tests for compatibility of genetic markers with the Hardy-Weinberg equilibrium (HWE) assumption as a means of quality control in genetic association studies (GAS) is to carry out this step of preliminary data analysis with the sample of non-diseased individuals only. We show that this strategy has no rational basis whenever the genotype-phenotype relation for a marker under consideration satisfies the assumption of codominance. Methods: The justification of this statement is obtained through rigorous analytical arguments. Results: The key result of the paper is that there holds the following proposition: under the codominance model, the genotype distribution of a diallelic marker is in HWE among the controls if and only if the same is true for the cases. Conclusion: The major practical consequence of that theoretical result is that under the codominance model, testing for HWE should be done both for cases and controls aiming to establish the combined (intersection) hypothesis of compatibility of both underlying genotype distributions with the HWE assumption. A particularly useful procedure serving this purpose is obtained through applying the confidence-interval inclusion rule derived by Wellek, Goddard, and Ziegler (Biom J. 2010; 52:253–270) to both samples separately and combining these two tests by means of the intersection-union principle.
1 Introduction
Testing the genotype frequencies observed for the SNPs involved for compatibility with Hardy-Weinberg equilibrium (HWE) condition is considered an important tool for the detection of technical issues and errors in genotyping. Accordingly, it has become one of the basic steps of quality control to be taken in the statistical analysis of any genetic association study (GAS), with genome-wide association study (GWAS) as a particularly important special case. Especially, in meta-analysis, one aims at controlling for deviations from HWE in the single studies to avoid systematic bias (see, e.g., Salanti et al. [1]). When focusing on case-control GAS, various recommendations for appropriate testing can be found in the literature. Since random mating is a mandatory assumption, most of them emphasize the necessity of HWE testing within the control group. A hybrid approach tailored for mixed populations comprising both cases and controls has been recommended by Wang and Shete [2]: These authors argue that testing for compatibility with HWE has in principle to be done for the total population which might include a substantial proportion of individuals being affected by the disease when the prevalence of the latter is large. In such a setting, testing for HWE with unaffected controls alone would exclude a relevant proportion of the total population, yielding potentially biased estimates of the proportions under evaluation. Wang and Shete [3] also address other potential sources of systematic deviations from HWE in the control group alone which are not considered technical issues. In GWAS, current recommendations focus on more stringent testing within the control group and suggest non-rigorous testing in cases (Turner et al. [4]; Marees et al. [5]). The tool EasyQC (Winkler et al. [6]) has become popular for genome-wide meta-analysis. It is to be used for filtering HWE violations in controls only, as suggested previously by Anderson et al. [7].
An issue with respect to which the practice of genetic epidemiology exhibits remarkable uniformity is the numerical coding of genotypes as determined in terms of single SNPs. In the vast majority of GAS, this is done through counting the allele hypothesized to be a potential risk factor for the disease under consideration (the allele of minor frequency in most cases). The corresponding genetic model, which is often called “codominance model,” reflects the notion that its impact on the (binary) phenotype is, on the usual log-odds scale, representable as a point on the straight line joining the points corresponding to the two homozygous genotypes. Thereby, it differs from the classical dominant and recessive models under which the effect of the heterozygous genotype equals that of one of the two homozygous genotypes. As is well known, the assumption of codominance is a presupposition for both of the most commonly used tests for association between a SNP and the disease under study (Cochran-Armitage test, logistic regression analysis with the allele count as the covariable of interest; for a comparison of both methods, see Wellek and Ziegler [8]).
The aim of this contribution is to show that the frequently recommended strategy of confining HWE testing to the controls has no rational basis whenever the genotype-phenotype relation for a marker under consideration satisfies the codominance model. This will be done through proving in a mathematically rigorous way the following logical equivalence relation: Under codominance, the distribution of a SNP-determined genotype is in HWE either in both underlying populations (i.e., cases and controls) or in none of them. In section 2, we start with introducing some notation and formulate the codominance model for an association study using cohort sampling. Afterward, we show how the model translates to the case-control setting which we are actually interested in. The main result is stated and proved in subsection 2.4 and section 3. In section 4, we outline and apply to some real examples a strategy of testing.
HWE compatibility in case-control studies on autosomal diallelic markers which is in line with the theoretical result stated in section 3 and also with the following crucial fact pointed out repeatedly in previous publications [9‒12]: A statistical testing procedure which is really suited for justifying the inclusion of a SNP in the final evaluation of the respective marker as a risk factor for the disease under study should have the form of an equivalence test tailored for establishing goodness rather than lack of fit of a marker with HWE. In the Discussion, it is pointed out that a qualitative difference between the results of the two tests for goodness-of-fit with HWE to be performed for each SNP in a GAS can be taken as a reason against basing the test for possible association of the SNP with the disease status on the customary codominance model.
2 Materials and Methods
2.1 Notation and Models
The data which are used in the statistical analysis of a genetic association study devoted to a binary phenotype and a diallelic marker (SNP) are the values observed for or previously assigned to the random variables appearing in a 2 × 3 contingency table of the generic form shown in Table 1. In order to derive the theoretical result making up the core of this contribution, it is convenient to represent the probabilities corresponding to the different cells of this contingency table as functions of the following population parameters:
Contingency table making up the frame of reference for the SNP-wise analysis of a genetic association study concerning a binary phenotype
. | Genotype . | |||
---|---|---|---|---|
Phenotype . | aa . | aA . | AA . | Σ . |
Control | X0 | X1 | X2 | N0 |
Case | Y0 | Y1 | Y2 | N1 |
Σ | S0 | S1 | S2 | N |
. | Genotype . | |||
---|---|---|---|---|
Phenotype . | aa . | aA . | AA . | Σ . |
Control | X0 | X1 | X2 | N0 |
Case | Y0 | Y1 | Y2 | N1 |
Σ | S0 | S1 | S2 | N |
τ = probability that an individual randomly chosen from the total population is affected by the disease under study (disease prevalence);
πdj = conditional probability of observing a genotype containing j A-alleles given the genotyped individual has been sampled from subpopulation d ∈ {0, 1}, with d = 0 and d = 1 indexing healthy and affected subjects, respectively, and j ranging over {0, 1, 2};
ξj = frequency of genotype j ∈ {0, 1, 2} in the total population.
The relationship between these parameters and the cell probabilities associated with the basic contingency table scheme displayed as Table 1 is made explicit in Table 2.
Cell probabilities associated with the contingency table shown as Table 1
. | Genotype . | |||
---|---|---|---|---|
Phenotype . | aa . | aA . | AA . | Σ . |
Control | π00 (1 − τ) | π01 (1 − τ) | π02 (1 − τ) | 1 − τ |
Case | π10τ | π11τ | π12τ | τ |
Σ | ξ0 | ξ1 | ξ2 | 1.00 |
. | Genotype . | |||
---|---|---|---|---|
Phenotype . | aa . | aA . | AA . | Σ . |
Control | π00 (1 − τ) | π01 (1 − τ) | π02 (1 − τ) | 1 − τ |
Case | π10τ | π11τ | π12τ | τ |
Σ | ξ0 | ξ1 | ξ2 | 1.00 |
2.2 Codominance Model for Studies Using Cohort Sampling
- 1.
The values s0, s1, s2 of the column totals S0, S1, S2 are considered as the sizes of three subsamples resulting from genotyping. The triplet (s0, s1, s2) not necessarily represents the genotype distribution in the total population.
- 2.
The Yj are mutually independent binomial random variables with parameters (sj, δj) where δj denotes the penetrance of genotype j, i.e., the disease risk of individuals carrying j A-alleles (j = 0, 1, 2).
- 3.
The codominance model assumes that, as measured on the log-odds scale, the increase in risk between “adjacent” genotypes remains constant. This is the case if and only if there exist constants α and β such that
2.3 Transforming the Codominance Model into a Model for Case-Control Studies
Equations (2) and (3) and the fact that the πdj are the same in the case-control study population as the total population justify the following conclusion: the codominance model, in its primary formulation referring to GAS using the cohort design, is satisfied if and only if the same holds true for the trinomial distributions to be analyzed in the corresponding case-control study.
2.4 Establishing Logical Equivalence between HWE Compatibility for Cases and Controls
with ηd: = πd1 /πd0, ζd: = πd2 /πd1.
The major conclusion from these expressions for the genotype frequencies π10 and π11 holding for the population of cases is as follows. If the parameter (π00, π01, π02) of the genotype distribution among the controls satisfy the HWE condition with risk-allele frequency p0 and the codominance model applies, then the population of cases is likewise in HWE with risk-allele frequency p1 = .
3 Results
The arguments and considerations made explicit in section 2 yield a rigorous proof of the fact that under the codominance model, in a genetic association study, the subpopulations to be compared with respect to the distribution of a diallelic marker are either both in HWE or in HW disequilibrium. Precisely speaking, the following statement turns out to hold true:
Remark, an alternative proof of the above result would be to derive it from Sasieni’s [13] Theorem 1.
4 Consequences for Preliminary HWE Testing in GeneticAssociation Studies
4.1 A Testing Strategy Taking the Result Derived in § 3 into Account
In the genetic epidemiology literature, the question of which of the two samples made available in a case-control association study should be checked by means of a preliminary statistical test for compatibility with HWE of the SNPs under evaluation has been a topic of controversial discussion (in addition to the references cited in the Introduction, see Wittke-Thompson et al. [14]; Ziegler [15], § 4.3.2). Essentially, the following two strategies are advocated:
- (I)
Before testing a diallelic marker for association with the phenotype, it should be made sure that it exhibits neither among the controls nor the cases a relevant amount of HW disequilibrium.
- (II)
Compatibility with HWE is a necessary requirement only for the population of controls. For cases, no test tailored for establishing the hypothesis of goodness-of-fit to HWE should be performed, since HW disequilibrium might indicate the existence of a true association being not just a spurious effect of population stratification or genotyping errors.
The result derived in this contribution admits the conclusion that strategy (II) is not really in line with the well-established practice of genetic epidemiology according to which genetic association studies are analyzed by fitting logistic regression models representing genotypes by single discrete covariates taking values 0, 1, or 2 rather than couples of binary dummy covariates. Our result implies that adopting this methodology and keeping it fully consistent requires that only SNPs which pass a preliminary test for goodness-of-fit with HWE both for controls and cases should be included in the confirmatory analysis of a GAS.
A recommendable, theoretically well-grounded strategy of testing the null hypothesis that a SNP exhibits relevant HW disequilibrium in at least one of the two subpopulations involved, against the alternative that the HWE condition is at least approximately satisfied both for cases and controls, proceeds in two steps.
Step 1: carrying out separately with each of the two samples a test at significance level α for goodness rather than lack of fit of the genotype distribution to the HWE model (for the theory and practical implementation of the exact version of this test, see [9, 10]. An asymptotic version based on confidence intervals which is particularly easy to implement, has been provided by Wellek, Goddard, and Ziegler [16]).
Step 2: acceptance of the SNP under consideration for evaluation as a potential genetic risk factor by means of logistic regression analysis if and only if each of the two individual tests performed in step 1 leads to rejecting the null hypothesis of relevant deviations from HWE.
4.2 Real Data Examples Illustrating the Use of the Proposed Strategy
termed relative excess heterozygosity by Wellek, Goddard, and Ziegler [16]. In this formula, the πj are generic symbols standing for the genotype frequencies in either of the two subpopulations, i.e., cases or controls. Following the proposal made in the same source, the equivalence margin for ω whose value under perfect HWE equals unity is set equal to 0.4, implying that the confidence interval has to be checked for inclusion in the acceptance interval (1/(1 + 0.4), 1 + 0.4) = (0.714, 1.400). In other words, the null hypothesis of relevant deviations from HWE has to be rejected if and only if it turns out that there holds both Cl (1 − α) > 0.714 and Cu (1 − α) < 1.400), where Cl (1 − α) and Cu (1 − α) denote the lower and upper confidence bound to ω at 1-sided confidence level 1 − α, respectively. As usual, the latter is chosen to be 95% throughout. All necessary computations are performed using the SAS/IML script ExactConfBound_RelHetZyg.sas to be found in [16] as accompanying material.
- 1.
An insertion/deletion polymorphism (INDEL) passing the preliminary goodness-of-fit test for both subpopulations.
The following data (absolute genotype frequencies; I: insertion, D: deletion) were obtained by Kee et al. [17] for a polymorphism of the angiotensin-converting enzyme (ACE) gene in an association study of risk factors for myocardial infarction:
Genotype . | II . | ID . | DD . | Σ . |
---|---|---|---|---|
X0j (Controls) | 193 | 292 | 215 | 801 (n0) |
X1j (Cases) | 184 | 460 | 217 | 861 (n1) |
Genotype . | II . | ID . | DD . | Σ . |
---|---|---|---|---|
X0j (Controls) | 193 | 292 | 215 | 801 (n0) |
X1j (Cases) | 184 | 460 | 217 | 861 (n1) |
Computing with these data the point estimates and asymptotic 95% confidence bounds for the parameter ω yields:
. | Cl (1 – α) . | . | Cu (1 – α) . |
---|---|---|---|
Controls | 0.8587 | 0.9646 | 1.0837 |
Cases | 1.0285 | 1.1510 | 1.2882 |
. | Cl (1 – α) . | . | Cu (1 – α) . |
---|---|---|---|
Controls | 0.8587 | 0.9646 | 1.0837 |
Cases | 1.0285 | 1.1510 | 1.2882 |
Obviously, both confidence intervals are contained in the interval (0.714, 1.400) so that the test for goodness-of-fit to the HW model rejects the null hypothesis of relevant deviations from HWE both for the controls and the cases. Thus, the INDEL under investigation satisfies the proposed criterion for being included in a subsequent test for association with the disease status.
- 2.
An SNP passing the preliminary goodness-of-fit test only for controls.
In an association study of genes predisposing to colorectal cancer, Shannon et al. [18] found the following absolute frequencies of genotypes of a SNP encoding in the methylenetetrahydrofolate reductase (MTHFR) enzyme:
Genotype . | TT . | CT . | CC . | Σ . |
---|---|---|---|---|
X0j (Controls) | 114 | 560 | 533 | 1,207 (n0) |
X1j (Cases) | 28 | 64 | 94 | 186 (n1) |
Genotype . | TT . | CT . | CC . | Σ . |
---|---|---|---|---|
X0j (Controls) | 114 | 560 | 533 | 1,207 (n0) |
X1j (Cases) | 28 | 64 | 94 | 186 (n1) |
With these data, the point estimates and asymptotic 95% confidence bounds for the parameter ω are computed to be:
. | Cl (1 – α) . | . | Cu (1 – α) . |
---|---|---|---|
Controls | 1.0179 | 1.1359 | 1.2676 |
Cases | 0.4755 | 0.6237 | 0.8182 |
. | Cl (1 – α) . | . | Cu (1 – α) . |
---|---|---|---|
Controls | 1.0179 | 1.1359 | 1.2676 |
Cases | 0.4755 | 0.6237 | 0.8182 |
The corresponding interval with endpoints Cl (1 − α) and Cu (1 − α) is included in the interval (1/1.40, 1.40) for controls but not for cases.
- 3.
A SNP passing the preliminary goodness-of-fit test only for cases.
The data for our next example are taken from a GAS of Alzheimer’s disease published by Bhojak et al. [19]. In the subgroup of APOE*4 carriers, these authors found the following absolute genotype frequencies of an IL-6 promoter polymorphism:
Genotype . | GG . | GC . | CC . | Σ . |
---|---|---|---|---|
X0j (Controls) | 34 | 22 | 12 | 68 (n0) |
X1j (Cases) | 110 | 148 | 44 | 302 (n1) |
Genotype . | GG . | GC . | CC . | Σ . |
---|---|---|---|---|
X0j (Controls) | 34 | 22 | 12 | 68 (n0) |
X1j (Cases) | 110 | 148 | 44 | 302 (n1) |
As point estimates and asymptotic 95% confidence bounds for ω, one obtains:
. | Cl (1 – α) . | . | Cu (1 – α) . |
---|---|---|---|
Controls | 0.3485 | 0.5446 | 0.8510 |
Cases | 0.8713 | 1.0637 | 1.2985 |
. | Cl (1 – α) . | . | Cu (1 – α) . |
---|---|---|---|
Controls | 0.3485 | 0.5446 | 0.8510 |
Cases | 0.8713 | 1.0637 | 1.2985 |
Clearly, only the second of the corresponding intervals (Cl (1 − α), Cu (1 − α)) satisfies the condition (Cl (1 − α), Cu (1 − α)) ⊆ (1/1.40, 1.40) so that approximate compatibility of the observed genotype distribution with HWE can be declared established at significance level α = 0.05 only for cases.
- 4.
A SNP passing the preliminary goodness-of-fit test for neither subpopulation.
Our final example is taken from a study of possible associations between polymorphisms of the cystatin C gene and Alzheimer’s disease in the Japanese population (Maruyama et al. [20]). The absolute frequencies found in this study for the CST3+73 G/A genotype were as follows:
Genotype . | GG . | GA . | AA . | Σ . |
---|---|---|---|---|
X0j (Controls) | 180 | 40 | 8 | 228 (n0) |
X1j (Cases) | 137 | 34 | 8 | 179 (n1) |
Genotype . | GG . | GA . | AA . | Σ . |
---|---|---|---|---|
X0j (Controls) | 180 | 40 | 8 | 228 (n0) |
X1j (Cases) | 137 | 34 | 8 | 179 (n1) |
As point estimates and asymptotic 95% confidence bounds for the parameter ω, one obtains with these data:
. | Cl (1 – α) . | . | Cu (1 – α) . |
---|---|---|---|
Controls | 0.3551 | 0.5270 | 0.7823 |
Cases | 0.3404 | 0.5135 | 0.7747 |
. | Cl (1 – α) . | . | Cu (1 – α) . |
---|---|---|---|
Controls | 0.3551 | 0.5270 | 0.7823 |
Cases | 0.3404 | 0.5135 | 0.7747 |
Since there holds neither (0.3551, 0.7823) ⊆ (1/1.40, 1.40) nor (0.3404, 0.7747) ⊆ (1/1.40, 1.40), the equivalence test rejects HW disequilibrium for neither subpopulation.
5 Discussion
In two out of the four examples considered in §4.2, the statistical analysis yields results which according to our basic theoretical result are inconsistent with the codominance model, due to the fact that the significance statements to be made in the test for goodness-of-fit with HWE for cases and controls are qualitatively different. Whenever this occurs, relying on the codominance assumption in evaluating the association of the genetic marker under consideration with the disease status can hardly be justified, and the recommendation to apply the intersection-union procedure described before loses its rational basis. The problem of how to combine the information on possible deviations from HWE obtained for both subpopulations has to be reconsidered, then. Developing a strategy of HWE testing for settings involving SNPs which do not satisfy the codominance model might be a worthwhile topic for future research.
Irrespective of whether or not the codominance model holds true, the combined testing procedure described in §4.1 maintains the predefined level α of significance as is guaranteed by the so-called principle of intersection-union (cf. R.L. Berger [21]; Wellek [11], §7.1). By this general result, it can be taken for granted that the risk of declaring some given SNP as evaluable with respect to its contribution to the risk of presenting the disease under investigation, although it fails to be in (approximate) HWE in at least one of both subpopulations, keeps below α (e.g., α = 0.05).
Researchers who do not see a problem in carrying out tests for lack of fit when establishing goodness of fit is obviously the aim, might be tempted to adopt the strategy outlined in section 4.1 with the traditional χ2 instead of the equivalence test as building block. This would not affect the justification for applying the recommended rule for combining the evidence in favor or against HWE contained in the data from the two arms of a GAS but would entail the risk of failing to keep the type-1 error risk of the combined testing procedure, at the target significance level. Actually, the risk of inferring approximate compatibility with HWE from a χ2-based p value being larger than the threshold α = 0.05 might grossly deviate from α to both directions. An instance where it is lower than required for making the decision based on the lack-of-fit procedure consistent with that taken in the goodness-of-fit test is the distribution obtained for the sample of cases of Example 4.2 (i). With X00 = 184, X01 = 460, X02 = 217, the traditional χ2 test for lack of fit yields a p value being as small as 0.0397, although the test for goodness-of-fit with the recommended equivalence interval (1/1.4, 1.4) decides in favor of approximate compatibility with HWE (at the same nominal significance level α = 0.05). Nevertheless, as a test for equivalence, the traditional χ2-based procedure will in general be anti-conservative, in the sense that the risk of declaring compatibility with HWE, although the true value of ω fails to fall between the limits 1/1.4 and 1.4, exceeds α. As has been shown by Wellek, Goddard, and Ziegler [16], using the equivalence interval (1/1.4, 1.4) leads to reasonable power to conclude compatibility with HWE with sample sizes of the order of magnitude available in standard GWAS.
Statement of Ethics
An ethics statement is not applicable because this study is based exclusively on analytical considerations and published literature.
Conflict of Interest Statement
The authors have no conflicts of interest to declare.
Funding Sources
There are no funding sources to acknowledge.
Author Contributions
S.W. formulated the problem addressed in this manuscript, developed the theory leading to its solution, worked out the proofs, and drafted an initial version of the manuscript. M.M.-N. and K.S. elaborated on the review of the existing literature on Hardy-Weinberg testing in genetic association studies and made many valuable suggestions for clarifying the conclusions to be drawn from the presented theoretical result.
Data Availability Statement
The data analyzed in this paper are publicly available and do not contain any information about individual subjects. The sources from which the data are taken are the references listed in the bibliography as items [17–20]. Further inquiries can be directed to the corresponding author.