key: cord-0280660-yhr7b7qy
authors: Zhang, Lin; Sun, Lei
title: Leveraging Hardy–Weinberg disequilibrium for association testing in case-control studies
date: 2020-11-16
journal: bioRxiv
DOI: 10.1101/2020.11.14.382796
sha: bfa0b9949631c7da47519d77cc4ceb6b2fcf05d9
doc_id: 280660
cord_uid: yhr7b7qy

In a case-control association study, deviation from Hardy-Weinberg equilibrium (HWE) or Hardy-Weinberg dis-equilibrium (HWD) in the control group is usually considered as evidence for potential genotyping error, and the corresponding SNP is then removed from the study. On the other hand, assuming HWE holds in the study population, a truly associated SNP is expected to be out of HWE in the case group. Efforts have been made in combining association tests with tests of HWE in the cases to increase the power of detecting disease susceptibility loci (Song and Elston (2006), Wang and Shete (2010)). However, these existing methods are ad-hoc and sensitive to model assumptions. Utilizing the recent robust allele-based (RA) regression model for conducting allelic association tests (Zhang and Sun (2020)), here we propose a joint RA test that naturally integrates association evidence from the traditional association test and a test that evaluates the difference in HWD between the case and control groups. The proposed test is robust to genotyping error, as well as to potential HWD in the population attributed to factors that are unrelated to phenotype-genotype association. We provide the asymptotic distribution of the proposed test statistic so that it is easy to implement, and we demonstrate the accuracy and efficiency of the test through extensive simulation studies and an application.

In a standard case-control study, deviation from Hardy-Weinberg equilibrium (HWE) or Hardy-28 Weinberg dis-equilibrium (HWD) in the control group is usually considered as evidence for po-29 tential genotyping error, and the corresponding SNP is then removed from the study. On the other sociation, and we demonstrate the accuracy and efficiency of the test through extensive simulation 52 and application studies. 53 2 Revisit the Robust Allele-based (RA) association test 54 The RA framework partitions a genotype as

if genotype is aa (0, 1) if genotype is Aa and c i = 0

(1, 0) if genotype is Aa and c i = 1

(1, 1) if genotype is AA

(1)

where c i iid ∼ Bernoulli(1/2). The term c i ensures the distributional assumption of the alleles; see 56 Section ?? for details.

Assume that we have an independent sample of size n. For a given biallelic SNP G and a 58 phenotype trait Y of interest, the RA regression framework for association analysis is

Based on (2), testing

evaluates the association between the SNP and phenotype of interest.

The RA regression (2) provides a unifying framework for allele-based association testing. In 62 the case of a binary outcome, the score test statistic of H 0 : β = 0 is T RA, binary, γ=0 = (p r −p s ) 2

( 1 2r + 1 2s )(p(1 −p) +δ )

,

whereδ =p 2 −p 2 is a sample estimate of the measure of HWD under H 0 . If we assume HWE 64 holds in the population (i.e. δ = 0), T RA, binary, γ=0, δ =0 is identical to the classical allelic test T allelic association. Recall the example that examines the association between HLA-DQ3 with cervical 84 intra-epithelial neoplasia (Sasieni, 1997) . The test of HWE yields a p-value of 3.54 × 10 −10 in the 85 controls and 0.0387 in the cases; see Table 1 for details. The HLA-DQ3 marker would have been 86 removed from the study following the standard HWE-based screening practice, but HLA-DQ3 is 87 well-known to be associated with cervical intra-epithelial neoplasia (Apple et al., 1994) . at the population level, with risk allele frequency p; see Table ? ? for definitions. The genotype 91 frequencies of (aa, Aa, AA) are then {(1 − p) 2 , 2(1 − p)p, p 2 }. In this work, we still use A to 92 denote the risk allele, but the risk allele can be either the minor allele with frequency ≤ 0.5 or 

and disease prevalence,

Thus, the genotype frequencies of AA and Aa in the case population are p (r) 2 = P(AA|Y = 1) = f 2 p 2 /K and p 1 /2 = f 2 p 2 /K + f 1 pq/K = ( f 2 p + f 1 q)p/K.

Thus, we expect departure from HWE in the case population, and the amount HWD is

Similarly, the amount of HWD in the control population is

Results in (7) and (8) additive genetic model ( f 1 = ( f 0 + f 2 )/2) and a recessive genetic model ( f 0 = f 1 < f 2 ).

Under the additive genetic model where f 1 = ( f 0 + f 2 )/2, (7) and (8) 

Thus, if the SNP is truly associated, i.e. f 2 = f 0 , δ ADD r and δ ADD s are both negative. The negative 115 sample estimates of δ r and δ s in Table 2 verify the theoretical derivations in (9). The test of HWE 116 is rejected in the case sample but not in the control sample at α = 0.05.

If the disease prevalence K is small, based on (9), the HWD in the control population, attributed 118 to true phenotype-genotype association, is too small to be detected using a control sample of mod-119 erate sample size at a stringent α level. Hence, if the genetic model is additive and the SNP is 120 in HWE at the population level, the HWD in the control population is expected to be negligible 121 without genotyping error. This may explain why testing for HWD in the control sample has been 122 used for detecting for potential genotyping error, as part of data quality control.

The second simulated example assumes a recessive genetic model where f 0 = f 1 < f 2 , with all 124 the other settings remain the same as the first example. Under a recessive genetic model, (7) and the two examples here also suggest that, rather than treating the HWD measure δ as a nuisance 136 parameter, a better method should utilize the expected HWD in the control sample (and the case 137 sample) attributed to true association. 138 We thus propose a unifying and flexible association test, RA joint test, that can be more pow-139 erful than the traditional association testing, by incorporating evidence for difference in HWD 140 between the case and control samples. Testing δ r = δ s instead of δ s = 0 (i.e. HWE holds in the 141 control sample) also provides robustness to departure from HWE in a sample due to genotyping 142 error or at the population level due to factors other than phenotype-genotype association. To develop the RA joint association association test for case-control studies, we adopt the general 145 RA regression framework in (2) but handle the δ parameter differently. For a case-control study 146 with independent samples, the RA joint association test is derived from

Unlike the original RA framework in (2) that assumes δ y being identical for all samples, the RA joint test allows δ y to be different for the case and control groups,

The RA joint association test evaluates 148 H 0 : β = 0 and δ r = δ s ,

and estimates δ r and δ s under the alternative hypothesis. We emphasize that in (11) we test if HWD 149 differs between the case and control samples, rather than if HWE holds in either the case or control 150 samples, to gain robustness to HWD due to genotype error or HWD at the population level.

The RA joint test statistic is

andβ =p r −p s . captures the differences in risk allele frequencies between the case and control samples; (ii) T δ 155 captures the difference in HWD estimates between the case and control samples; (iii) T p,δ adjusts 156 for the correlation between T p and T δ ,

.

Note that T p is similar to T allelic, Schaid in (??) (Schaid and Jacobsen, 1999), except that T p subtracts 161 cov(p r −p s ,δ r −δ s ) 2 / var(δ r −δ s ) in the denominator to account for the dependence between T p 162 and T δ .

Assume the case and control samples are independent of each other as in a standard case-control 164 study, it is easy to show that var(δ r −δ s ) = var(δ r ) + var(δ s ), cov(p r −p s ,δ r −δ s ) = cov(p r ,δ r ) + 165 cov(p s ,δ s ), and var(p r −p s ) = var(p r ) + var(p s ), where the variance-covariance components can 166 be estimated in the case and control samples separately.

For the case sample, as the sample size r → ∞, the asymptotic distribution of (p r ,δ r ) T is

Σ s is defined similarly by replacing r in (14) with s; see Supplementary Material ?? for the detailed 170 derivation.

The variance-covariance matrix of (β ,δ r −δ s ) T can then be written as The traditional 1 d.f. additive test codes the genotypes aa, Aa and AA as 0, 1 and 2, respectively.

The additive test for a case-control study uses logistic regression,

and evaluates the association between Y and G by testing

The traditional 2 d.f. genotypic test treats the three genotypes aa, Aa and AA as categories. The

genotypic test for a case-control study also uses logistic regression,

where, without loss of generality, genotype aa is chosen to be the baseline category, evaluates the association between Y and G.

For fair comparison with the proposed RA joint test, we evaluate both (17) and (19) Suppose the risk allele frequency of G is p with population level departure from HWE, δ = p 2 − p 2 .

The population level genotype frequencies of aa, Aa and AA are

We follow the same sampling procedure used for the two illustrative examples in Section 3, 208 but without the assumption of HWE at the population level. We sample the genotypes based on 209 probabilities {p 0 , p 1 , p 2 }, and each genotype G = k is assigned to the case group with probability Table 3 summarizes the empirical type I error rate of the proposed RA joint test for each 219 combination of p and δ , which shows that the proposed test is accurate. Figure 1 Table 4 : Six common genetic models and their corresponding penetrance specifications for the power studies. For additive, additive on log-scale, recessive, and dominant genetic models, we assume the homozygote relative risk f 2 / f 0 = 1.5 and the corresponding f 1 is determined by the genetic model assumptions; for heterozygous advantage genetic model, we assume the heterozygote relative risk f 1 / f 0 = 1/1.5, and f 0 = f 2 ; for heterozygous disadvantage genetic model, we assume the heterozygote relative risk f 1 / f 0 = 1.5, and f 0 = f 2 .

can be capped at 11% based on a theoretical comparison between χ 2 1 and χ 2 2 distributions, re- and control allele frequencies using an adjusted variance.

The power study in Figure 2 assumes HWE holds in the population. We show in Supplemen-

tary Figure ? ? that, as the population level HWD (δ ) increases, the power of the traditional 1 d. we follow the same simulation setting as in Section 6.2, except that we use α = 5 × 10 −8 and 268 increase the sample size to r = s =2,500 for both the case and control groups. The empirical type 269 1 error is evaluated based on 10 10 replicates for each parameter setting. well-controlled type 1 error rate and is robust to departure from HWE at the population level. When 272 the risk allele frequency is very low or very high (e.g. 0.1 or 0.9), the test is slightly conservative.

On the other hand, the traditional 2 d.f. genotypic test appears to be generally conservative at the 274 α = 5 × 10 −8 level, which has not been noted before to the best of our knowledge. 275 We then assess the power of the proposed RA joint test and the traditional 2 d.f. genotypic test. 276 We follow the same simulation setting as in Section 6.3 using the six common genetic models; 277 see Table 4 for details. We fix f 2 = 0.1 and vary the heterozygote relative risk or the homozygote contains a sample of 2,500 cases and 2,500 controls. However, it appears that the power of T RA joint is practically identical to that of T genotypic .

In the next section, we provide theoretical insights on the power of the two 2 d.f. association tests In this section, we examine the theoretical power of the two tests through non-central parameters. 301 We first show that the general form of the non-central parameter of T genotypic is 302 NCP genotypic = r

where

and π 01 = s/r is the ratio of the number of controls to the number of cases. Due to the complex has shown that the power difference between a χ 2 1 and a χ 2 2 test statistic can be capped at 11%. We Tables 6 and ? ?.

In the empirical power studies, it appears that the proposed RA joint test is practically identical 360 to the traditional genotypic test. We thus further compare the theoretical power between the two 361 tests (both follow the χ 2 2 distribution) via non-central parameters. We show that the non-central 362 parameter of the proposed RA joint test is larger than that of the traditional genotypic test under 363 all six genetic models in balanced case-control studies. When the effect size is small, the dif-364 ference between the NCPs of the two tests are marginal, but NCP RA joint is slightly higher than 365 NCP genotypic ; see Supplementary Figure ? ?. As the effect size increases, the difference in NCPs 366 between the two tests can be substantial, indicating that T RA joint can be more powerful to detect 367 truly associated SNPs in GWAS, which is also verified by the application studies.

In the simulated and theoretical power studies, we have assumed equal numbers of cases and 369 controls, i.e. π 01 = s/r = 1. However, in real-life genetic association studies, often there are more 370 controls than cases (i.e. π 01 > 1). We show that the effect of π 01 on NCP RA joint and NCP genotypic 371 depends on the genetic models, risk allele frequencies and the disease penetrance probabilities; see 372 Supplementary Figure ? ? for details.

When the frequency of the risk allele is close to 0 or 1, e.g. p A < 0.05 or p A > 0. 95 By convention, allele A is usually referred as the minor allele (i.e. p A ≤ 0.5) with the assump-386 tion that diseases are rare, and it is the minor allele that increases the risk. In this work, we allow 387 p A to be bigger than 0.5, and broaden the notion of 'disease' traits to any rare or common phe-388 notype traits. For example, in the current genetic studies of susceptibility and disease severity of 389 COVID-19, the risk allele can be a major allele for disease susceptibility which is common; the 390 risk allele is likely to be the minor allele for disease severity because it is rare (Lin et al., 2020).

Therefore, in the simulated and theoretical power studies of this work, the range of p A is from 0.05 392 to 0.95.

From this whole range of p, we observe the power of the proposed RA joint test, the traditional

Data quality control in genetic case-control association studies

Hla dr-410 dq associations with cervical carcinoma show papillomavirus-type specificity

The uk biobank resource with deep phenotyping and genomic 413 data

The x factor: A robust and pow-415 erful approach to x-chromosome-inclusive whole-genome association studies

Estimation of significance thresholds for genomewide 420 association scans. Genetic Epidemiology: The Official Publication of the International Genetic 421

On statistical power for case-control host genomic studies of 423 covid-19

A tutorial on conducting genome-wide association studies: Quality control and statistical 426 analysis

From genotypes to genes: doubling the sample size

Biased tests of association: comparisons of allele fre-429 quencies when departing from hardy-weinberg proportions

A powerful method of combining measures of association and 432 hardy-weinberg disequilibrium for fine-mapping in case-control studies

A tail strength measure for assessing the overall univariate 435 significance in a dataset

A test for genetic association that incorporates information about 437 deviation from hardy-weinberg proportions in cases

Rational inferences about departures 440 from hardy-weinberg equilibrium