|Home | About | Journals | Submit | Contact Us | Français|
Testing Hardy-Weinberg Equilibrium (HWE) in the control group is commonly used to detect genotyping errors in genetic association studies. We propose a likelihood ratio test for testing HWE in the study population using both case and control samples. This test incorporates underlying association models. Another feature is that, when we infer the disease-genotype association, we explicitly incorporate HWE or a possible departure from Hardy-Weinberg Equilibrium (DHWE) into the model. Our unified framework enables us to infer the disease-genotype association when a detected DHWE needs to be part of the model after causes for the DHWE are explored. Real datasets are used to illustrate the application of the methodology and its implication in genetic association studies. Our analysis and interpretation touch on genotyping errors, population selection, population stratification, or the study sampling plan, all delicate issues that could be the cause of DHWE.
Hardy-Weinberg Equilibrium is used to describe the genotype distribution of a population when it is large, self-contained, and randomly mating. The equilibrium can be summarized as, if p is the frequency of one allele (A) and q is the frequency of the alternative allele (a) for a biallelic locus, then the HWE-expected frequency will be p2 for the AA genotype, 2pq for the Aa genotype, and q2 for the aa genotype. The three genotypic proportions should sum to 1, as should the allele frequencies (Hardy, 1908; Weinberg, 1908).
Many methods have been developed to test HWE. Weir (1996) and Emigh (1980) provide summaries of these methods. A χ2 test is commonly used to assess a departure from HWE. Exact tests of a departure from HWE have been developed for studies with small sample sizes (Haldane 1954; Wigginton, Cutler, and Abecasis, 2005).
Testing HWE is commonly conducted for genotyping quality control (Gomes, Collins, et al., 1999; Xu, Turner, et al., 2002). Some view the testing as an essential step in genetic association studies (Xu, Turner, et al., 2002; Thakkinstian, McElduff, et al., 2005); however, others caution such use (Nielsen, Ehm, and Weir, 1999; Wittke-Thompson, Pluzhnikov, Cox, 2005; Zou and Donner, 2006). Nielson, Ehm, et al. (1999) point out that HWE is generally expected to be distorted in the case sample in the region of association. Zou and Donner (2006) suggest testing for HWE should not be used as a tool for identifying genotyping errors when it is tested in a single sample. Wittke-Thompson, Pluzhnikov, Cox (2005) provide a framework to guide the interpretation of a DHWE for case-control studies. They suggest that if a DHWE in cases or in both cases and controls is detected, it does not necessarily imply genotyping errors. Rather than discarding the data, the underlying disease-genotype association should be investigated. The association may explain the observed DHWE. If not, other possible explanations such as “genotyping error, chance, failure of assumptions underlying Hardy-Weinberg expectations” should be explored. In their framework, they explicitly assume the genotype is in HWE in the population.
The current work develops a likelihood ratio test for testing HWE in the study population. Our test uses data from both case and control samples, and the procedure accounts for the underlying disease models. We estimate the parameters in the model by minimizing the deviance, thus, the estimates are maximum likelihood estimates. The difference between the deviances of two nested models follows a χ2 distribution. This forms a likelihood ratio test for the population HWE.
When we infer the disease-genotype association, HWE or a possible DHWE are explicitly incorporated into the model. A DHWE could be due to a variety of reasons such as genotyping errors, population selection, population stratification, or the sampling plan of the study. If they come into play, then it is likely these problems will impact both cases and controls. The purpose of the likelihood ratio test of the population HWE is to prompt the investigators to think about these issues in addition to possible genotyping errors. If a DHWE is detected and genotyping errors are ruled out, our unified framework enables us to model the association under DHWE in contrast to the current practice that testing HWE is a middle step in the analysis of association studies.
The rest of the manuscript is organized as follows. Section 2.1 summarizes common genetic disease models. In Section 2.2, we develop the likelihood ratio test for population HWE. In Sections 3.1 and 3.2, we demonstrate our methods in detail using data from two genetic association studies conducted at Vanderbilt University Medical Center. Then we revisit some examples discussed by Wittke-Thompson, Pluzhnikov and Cox (2005). We developed the analysis software in R and it can be obtained from the corresponding author.
We first summarize the common disease models presented in Wittke-Thompson, Pluzhnikov, Cox (2005). Then we develop a likelihood ratio test for testing HWE in the study population.
Wittke-Thompson, Pluzhnikov, Cox (2005) explicitly assume the susceptibility locus is in HWE in the study population. This assumption implies the genotype distribution follows
where p is the population frequency of the wild-type allele (A), and q is the population frequency of the disease-susceptibility allele (a).
Let α be the baseline disease penetrance in homozygotes (AA), and Y = 1, 0 denote the study outcome diseased or not. They present the following general disease model:
where β is the relative risk of disease for the hetrozygotes Aa in reference to homozygotes AA and γ is the relative risk of disease for the homozygotes aa.
The prevalence of disease in the population is restricted as
It is recognized that the disease prevalence Kp can not be estimated in a case-control study, thus it has to be obtained through external studies.
In this framework, Wittke-Thompson, Pluzhnikov, Cox (2005) propose to estimate the parameters θ = (β, γ, q) by minimizing the general χ2 statistic. α is obtained through constraint (3) after the estimates (, and ) are obtained.
In this section, we first expand the common disease models to incorporate whether or not the genotype is in HWE in the population. Then we develop a test for the population HWE expressed as a null hypothesis that the susceptibility locus is in HWE in the population versus the alternative hypothesis that the susceptibility locus is not in HWE in the population. Here we define the study population as the population from which cases and controls are drawn and to which study findings will be generalized.
The null hypothesis H0 is expressible as (1), and under H0 the disease models are described in Section 2.1. The alternative hypothesis Ha can be expressed as the genotype distribution in the population
where p2 = 1 − p0 − p1. Under Ha, Kp is constrained as
Table 1 summarizes a typical data set from a traditional unmatched case-control study. The objective is to explore the relationship between disease status and a single-locus 2-allele genotype denoted as AA, Aa, and aa.
Assuming common disease model (2), under either H0 or Ha, the conditional distributions of the genotype in cases and in controls are
Conditional on the study design parameters, n.1 cases and n.2 controls, the expected numbers of genotypes in cases are:
The expected number of genotypes in controls can be expressed as replacing n.1 with n.2 and the probability expressions with (7).
We fit models by minimizing the deviance function
over the parameter space of θ, where θ depends on the specific model. Due to the limitation that the data has only 4 independent observations, only models more parsimonious than the general models can be estimated under the alternative hypothesis. General models under Ha are saturated with θ = (p0, p1, β, γ). In Section 3, we will only consider examples for which recessive disease models (i.e. β = 1) or dominant disease models (i.e. β = γ) are appropriate.
We use the R function nlm to obtain parameter estimates in the models. The standard error of the estimates is obtained through the inverse of the Hessian matrix and through the delta method since proper transformations of the parameters are needed to facilitate the algorithm.
Under H0 and Ha, we fit the same disease model (2), and the difference of the deviances follows a χ2 distribution with one degree of freedom under regularity conditions (Section 4.5.4 of Agresti, 2002, page 141−142). This forms a likelihood ratio test for the population HWE, H0 versus Ha.
In this section, we illustrate the application of our methodology and its implication in genetic association studies using several real datasets. The example in Section 3.1 explores a set of possible explanations for the observed DHWE. Section 3.2 discusses the potential impact of the sampling plan on the analysis. For both examples, our analysis results always prompted the investigators to reevaluate genotyping quality. After genotyping errors were ruled out, additional issues were explored as possible explanations for the DHWE in the study population. Consequently, the inference of the disease-genotype association would be made under the best-fit disease model.
We applied our methodology to an association study involving genetic variations in an accessory chloride channel subunit and hypertension in the Ghanaian population. BSND-V43I was identified as a common polymorphism in the non-Caucasian population. Functional examination of this variant demonstrated a partial loss-of-function variant when heterologously expressed with ClC-Kb in cultured cells. The BSND-V43I genotypes (GG, AG, AA) are (155, 27, 4) in the cases and (408, 55, 15) in the controls (Sile, Gillani, et al., (2007)). In current practice, where HWE is separately tested in cases and in controls, both demonstrate a significant departure from HWE (p = 0.043, and p < 0.001, respectively). Applying the general disease model, the best fit is a model with α = 0.2006, β = 0.99, γ = 0.85, and q = 0.11 using the essential hypertension prevalence of Kp = 0.20 (based on World Health Organization report (2005)). Here GG is the reference group and q is the A allele frequency in the population. The best-fit model has a lack-of-fit statistic of 36.66 with 1 degree of freedom, exhibiting a significant lack of fit.
The estimated β indicates a recessive disease model. Therefore, we separately fit a recessive model under H0 and Ha. Results are presented in Table 2. The recessive model under H0 had a similar fit to the general model. A note should be made that the deviance for the general model under H0 is actually the general χ2 statistic and it is different from the deviance defined in (8) when we fit our models. We did that since we would like to analyze this data set using both the methods of Wittke-Thompson, Pluzhnikov and Cox (2005) and our methods. This explains why the deviance for the recessive model under H0 is smaller than the general model under H0. The deviance indicates the recessive model under Ha offers significant improvement over the same model under H0. The recessive model under Ha no longer demonstrates a significant lack-of-fit. The likelihood ratio test suggests significant evidence against H0. Hence, we conclude that the genotype is not in HWE in the Ghanaian population from which the cases and the controls were drawn.
There are several possible explanations for the departure from HWE in this example with the BSND-V43I allele. Some of the explanations involve genotyping errors, population stratification, sampling methods, selection, non-random mating, and possibly chance.
Historically, deviation from HWE has been attributed to genotyping errors; however the investigators had several measures in place to avoid such problems. These measures included unstructured sample-numbering system regarding cases and controls, blanks and sequence-verified controls in each plate (Sile, Gillani, et al., 2007).
In addition, deviation from HWE could stem from population stratification or admixture. However, previous studies by Adeyemo, Chen, et al., (2005) demonstrated that these issues are negligible in the Ghanaian population that Sile, Gillani, et al., (2007) examined.
Furthermore, sampling methods could not explain this observation of deviation from HWE. The samples were not ascertained with regards to certain clinical or genetic phenotype or prior knowledge of any disease status. The exclusion criterion was the presence of an acute illness. We believe the cases and the controls are random samples.
Finally, selection appears to be a possible explanation. Unfortunately, there is no data regarding other markers in linkage disequilibrium (LD) in this region to further examine the issue of selection. Using electrophysiology patch clamping, Sile, Gillani, et al., (2007) demonstrated that BSND-V43I is a partial loss-of-function polymorphism. They concluded, based on their functional data, that susceptible subjects with this allele might be protected from developing hypertension. Based on their conclusion we think that selection might be an explanation for the observed deviation from HWE. Further studies examining patterns of variations and LD are needed to determine if selection is indeed present in this region. Additional possible explanations include non-random mating and chance, for which we can not evaluate in this study.
We also applied our methodology to an association study of TGFβ1 Codon 10 Polymorphism and Familial Pulmonary Arterial Hypertension (FPAH). Dr. John Phillips III and his group genotyped TGFβ1 Codon 10 on a cohort of 120 FPAH patients (probands) and 51 of their relatives. Every person in the cohort has the BMPR2 mutation. The hypothesis is that codon 10 T to C transition increases expression and circulating levels of TGFβ1, thus increasing the chance for FPAH. The TGFβ1 codon 10 SNP genotypes (TT, CT, CC) are (29, 78, 13) in the cases and (17, 28, 6) in the controls (Phillips III, Poling, et al., 2008). When tested separately, cases demonstrated a significant departure from HWE (p = 0.0004), but controls did not. Fitting the general disease model, the best fit is a model with α = 0.0000066, β = 2.06, γ = 1.05, and q = 0.39 using a FPAH prevalence of Kp = 0.00001 (Online Mendelian Inheritance in Man, 2007). Here TT is the reference group, β and γ are the relative risk of CT and CC groups, respectively, and q is the C allele frequency in the population. This best-fit model has a lack-of-fit statistic of 1.17 with 1 degree of freedom, indicating no significant lack of fit. We conclude that the underlying disease model explains the observed departure from HWE in the cases, though it is difficult to interpret the model with β = 2.06 and γ = 1.05.
Based on preliminary analyses conducted by Dr. John Phillips III and his group, it appears that the association between this genotype and FPAH follows a dominant disease model. We separately fit a dominant model under H0 (i.e. in the population the TGFβ1 Codon 10 SNP genotypes are in HWE) and Ha (i.e. in the population the TGFβ1 Codon 10 SNP genotypes are not in HWE). The results are summarized in Table 3. The dominant model does not indicate significant lack-of-fit under either H0 or Ha. The likelihood ratio test, however, suggests significant evidence against H0. Hence, we conclude that the genotype is not in HWE in the population from which the cases and the controls were drawn. We communicated this finding to the investigator. The investigation team first regenotyped this SNP; and genotyping errors were ruled out. One possible explanation for the DHWE in the population might be the plan of the study ascertainment. FPAH is a rare disease, and it is established that people with the BMPR2 mutation are at a higher risk for this disease. In the study sample, every subject has the BMPR2 mutation. Thus, the study sample comprises a stratum from the general population, rather than a representative sample. Secondly, the controls are relatives of the probands. This works against the principles of case-control studies. In this example, however, the primary objective, given the well-established risk factor of the BMPR2 mutation, is to evaluate the additional effect of TGFβ1 Codon 10 SNP on this rare disease. Thus the design is appropriate. We recognize that since they are related, the controls might have the tendency to show a similar distribution of TGFβ1 Codon 10 genotype to the diseased. However, any bias introduced would be toward the null and makes the investigators’ findings conservative.
Among the examples discussed by Wittke-Thompson, Pluzhnikov and Cox (2005), we focus on the studies that showed significant lack-of-fit by their best-fit recessive models. The study by Ozaki, Ohnishi, et al. (2002) had the heterozygote relative risk estimate of 1.024, so it is also included. Under the assumption that the genotype is in HWE in the population, the lack-of-fit by the best-fit genetic disease models suggests the genotype disease association is an unlikely explanation for the observed DHWE in patients or in controls.
We first fit recessive models, i.e., with β fixed at 1, under the hypothesis that the genotype is in HWE for these examples. Table 4 summarizes the results. The parameter estimates and the goodness-of-fit statistics are almost identical or very close to those reported in Table 1 of Wittke-Thompson, Pluzhnikov and Cox (2005). All the best-fit recessive models show significant lack-of-fit.
We then fit the same recessive models, but under the assumption that the population is not in HWE. Results are reported in Table 5. With two parameters, p0 and p1, to represent the genotype frequencies, the allele frequency q is calculated as p1/2 + (1 − p0 − p1). The estimates of q and α are similar to those in Table 4, and some estimates of γ have changed significantly. As expected, the deviances are now smaller, and all but one demonstrate significant improvement over the models reported in Table 4. For the study with PubMed ID 11889073 (see Table 5), the recessive model still shows significant lack-of-fit and the investigator needs to look for possible reasons other than the DHWE in the population for an explanation of the lack-of-fit.
Comparing the models under H0 and Ha using the likelihood ratio test, we see that much of the lack-of-fit shown in the recessive disease models (Table 4) can be explained by assuming the genotypes are not in HWE in the population. The likelihood ratio tests are reported in Table 6.
These findings raise questions for the investigators about the reason for the departure from HWE in their study population. Sections 3.1 and 3.2 present two detailed examples in which a set of possible explanations were explored to explain the departure. In addition to genotyping errors, we suggest the investigators look into similar issues for possible explanations of the DHWE. On the other hand, whether the genotype is in HWE in the study population plays an important role in making inference about the genotype disease association. Therefore, this assumption should be assessed explicitly in the model.
HWE is commonly tested separately in cases and in controls for genotyping quality control. Several researchers have expressed concern about this practice (Nielsen, Ehm, and Weir, 1999; Wittke-Thompson, Pluzhnikov, Cox, 2005; Zou and Donner, 2006). The current work proposes a likelihood ratio test for testing HWE using both case and control samples. If the problems that HWE testing is intended to detect, such as genotyping errors or population stratification, come into play, these problems are likely to impact both cases and controls. Rather than the current approach of separately testing HWE in cases and in controls as a middle step in the analysis of association studies, our methods test HWE in the study population. We explicitly incorporate HWE or a possible DHWE into the model when we infer the underlying disease-genotype association. If genotyping errors are ruled out and the DHWE is plausible, our methods provide a means to study the association. The observation that some of the estimates of γ in Table 5 changed significantly from the estimates in Table 4 underlines the message that the association estimates depend on the assumption.
Testing HWE in the study population also has implications beyond genotyping quality control. Some genetic methods depend on the assumption that the population is in HWE (Cheng and Chen (2005) and Cheng and Lin (2005)). Our methods provide investigators a useful tool to test HWE in the study population before they apply the analysis methods in those settings.
The examples in Section 3 illustrate the application of the methodology and its implication in genetic association studies. As demonstrated in the examples, a DHWE in the study population could be due to reasons other than genotyping errors, such as population stratification, population selection, the study sampling plan, or failure of the assumptions underlying HWE. Although our methods appear to carry a simple message, they touch on these delicate issues. We suggest a detected DHWE in the study population be investigated with this in mind.
One limitation of our methods is that we assume the disease prevalence is fixed. It is recognized the disease prevalence can not be estimated in a case-control study, thus it has to be obtained from external sources. In the example of Section 3.1, we also analyzed the data using essential hypertension prevalence of K = 0.15 and K = 0.25. Results were similar to the reported and therefore the conclusion is unchanged. In practice, we suggest investigators conduct sensitivity analyses using a range of plausible estimates of disease prevalence.
We acknowledge John Phillips III, Scott Williams, Chun Li, and Daniel Zelterman for helpful discussions. We also acknowledge John Phillips III for giving us access to the data set of TGFβ1 codon 10 polymorphism and familial pulmonary arterial hypertension association study. This research was supported in part by the U.S. National Institutes of Health with grant RR00095 awarded to the General Clinical Research Center at Vanderbilt University Medical Center.