Robust ranks of true associations in genome-wide case-control association studies

Abstract

In whole-genome association studies, at the first stage, all markers are tested for association and their test statistics or p-values are ranked. At the second stage, some most significant markers are further analyzed by more powerful statistical methods. This helps reduce the number of hypotheses to be corrected for in multiple testing. Ranks of true associations in genome-wide scans using a single test statistic have been studied. In a case-control design for association, the trend test has been proposed. However, three different trend tests, optimal for the recessive, additive, and dominant models, respectively, are available for each marker. Because the true genetic model is unknown, we rank markers based on multiple test statistics or test statistics robust to model mis-specification. We studied this problem with application to Problem 3 of Genetic Analysis Workshop 15. An independent simulation study was also conducted to further evaluate the proposed procedure.

Background

For a large genetic study, a two-stage analysis is often employed. At the first stage, each marker is tested for association with a disease. The p-values of all markers are ranked. Then some of the most significant markers are analyzed in the second stage. This two-stage analysis reduces the number of hypotheses to be tested in the second stage. Hence, it enhances the power to identify true marker susceptibility to the disease. However, it is important to know how many of the most significant markers one should study in the second stage so that the probability that one or several true markers will be studied in the second stage is greater than a given value. On the other hand, when a given number of the most significant markers is selected, it is important to know the probability that this list of markers would contain one or more true markers. A small list of the most significant markers may not contain any true markers at all, which leads to spurious associations or negative findings in the second stage.

Zaykin and Zhivotovsky [1] used p-values of a single test statistic to rank markers. In a case-control study for complex diseases, three trend tests can be applied under the recessive, additive, and dominant models. Because the genetic model of the marker is uncertain, ranking the markers with a single test statistic may not be robust when another genetic model is correct. Using the first simulated data set of Problem 3 from Genetic Analysis Workshop (GAW) 15, we study robust ranking when the underlying genetic model is unknown and examine whether robust test statistics would lead to robust rankings of about 10 K single-nucleotide polymorphisms (SNPs). The properties of the proposed robust ranking procedures are then further examined by an independent simulation study.

Methods

Notation and model

Consider a SNP with alleles D and d and frequencies p and q = 1 - p, respectively. In a case-control design, r cases and s controls are independently sampled from a population. The genotype counts of three genotypes G₀ = dd, G₁ = Dd, and G₂ = DD are denoted as (r₀, r₁, r₂) in cases and (s₀, s₁, s₂) in controls, which follow multinomial distributions mul(r: p₀, p₁, p₂) and mul(s: q₀, q₁, q₂), respectively. Denote the disease prevalence as k and penetrances as f _i = P(case|G _i ) for i = 0, 1, 2. By the Bayes Theorem, p _i = g _i f _i /k and q _i = g _i (1 - f _i )/(1 - k), where g _i = P(G _i ). Without loss of generality, assume that D has high risk. Then the null hypothesis of no association can be stated as H₀: f₀ = f₁ = f₂ = k. The alternative hypothesis is H₁: f₀ ≤ f₁ ≤ f₂ with at least one inequality. The genotype relative risks (GRRs) are defined as λ₁ = f₁/f₀ and λ₂ = f₂/f₀. The recessive, additive, and dominant models are referred to as λ₁ = 1, λ₁ = (1 + λ₂)/2, and λ₁ = λ₂, respectively [2–4].

Trend tests and robust tests

To test association using case-control data, the Cochran-Armitage trend test (CATT) has been proposed [2–4], which can be written as

Z_x=\frac{n^{1/2}{\displaystyle \sum _{i=0}^2{x}_i\left(s{r}_i-r{s}_i\right)}}{{\left[rsn\left\{n{\displaystyle \sum _{i=0}^2{x}_i^2{n}_i-{\left({\displaystyle \sum _{i=0}^2{x}_i{n}_i}\right)}^2}\right\}\right]}^{1/2}},

(1)

where (x₀, x₁, x₂) = (0, x, 1) and 0 ≤ x ≤ 1. Given x, Z _x follows asymptotically N(0,1). The choice of x is 0, 1/2, and 1 for the recessive, additive/multiplicative, and dominant models, respectively [5]. In practice, however, the true genetic model is unknown. Hence the robust tests, maximin efficiency robust test (MERT) and maximum test (MAX), can be applied, which are given by MERT = (Z₀ + Z₁)/{2(1 + ρ)}^1/2 and MAX = max(|Z₀|, |Z_1/2|, |Z₁|), where ρ = [n₀n₂/{(n₀ + n₁)(n₁ + n₂)}]^1/2 [4]. Note that Pearson's association test can also be used. However, Zheng et al. [6] showed that the MAX is often more powerful than the Pearson chi-squared test for a case-control design. Comparison of MERT and MAX can be found in Freidlin et al. [7]. The MAX and MERT have also been applied to other designs for GAW14 [8, 9].

Ranking markers with multiple statistics

When the genetic model is unknown, the three CATTs (Z₀, Z_1/2, Z₂) are calculated for each of M SNPs. Then the p-values of MERT and MAX can be obtained for ranking. However, computing the p-value of MAX needs extensive simulation. Thus, alternatively, the minimum of the p-values (min p) of the three CATTs can be used for ranking. Rather than ranking M SNPs based on any single CATT, we propose ranking the SNPs by the MERT and the minimum of the p-values. We expect that ranking SNPs based on this approach would be more robust compared to ranking by a single CATT when the ranks by the three CATTs are quite different.

Results

Application to GAW15

As an application, we consider the first simulated data set of Problem 3 from GAW15. A simulated data set was considered, as we knew that there were eight candidate genes. One of them at chromosome 6 with physical location 32,484,648 bp was simulated based on the DRB1 locus of the HLA gene. We selected four SNPs closest in physical distance to the eight known candidate genes as candidate SNPs. We examined the ranks of the 32 candidate SNPs among all 9187 SNPs. All 2000 unrelated controls were used. For the affected sib-pair (ASP) data, we selected an affected sib (case) with the first individual ID from each family. A total of 1500 unrelated cases were used. In the simulated data set, genotypes of all 9187 SNPs from 22 chromosomes were generated (no missing genotypes and no genotyping errors). All SNPs had minor allele frequency (MAF) greater than 1% and there were no monomorphisms. Because we considered the CATTs, Hardy-Weinberg equilibrium in the population was not required [2]. If any genotype count in cases or controls was 0, 0.5 was added to all genotype counts in cases and controls.

After Bonferroni correction for Z₀ (Z_1/2 and Z₁), there were 5 (7 and 7) SNPs among the 32 candidate SNPs that had Bonferroni-corrected p-values less than 0.05. All three CATTs, the MERT, and the minimum of the p-values of the three CATTs were used to rank all 9187 SNPs. The ranks of the 32 candidate SNPs are reported in Table 1 by five different ranking methods. The results are summarized below: 1) in the candidate gene DRB1 of HLA (chromosome 6, location = 32,484,648), four of the six most significant candidate SNPs are in this region. This implies that when the sample size and a genetic effect are large, a strong candidate gene should contain several SNPs at the top of the list of most significant SNPs. 2) Using a single CATT to rank SNPs may not be robust, and using MERT or the minimum p-value is more robust. For example, the SNP (chromosome 6, location = 37,363,880) has rank of 6 using either Z_1/2 or Z₁, and 8172 when Z₀ is used. But the ranks of this SNP by MERT and minimum p-value are 10 and 6, respectively. 3) When the ranks by the three CATTs are quite different, the ranks by the robust methods are usually in the middle. 4) With a sample size of 3500, some candidate SNPs have ranks larger than those of null SNPs. Thus, selecting only the most significant SNPs from the genome-wide scan for further analysis may exclude some true associations or candidate genes. This information is particularly important for cost-efficient two-stage design for genome-wide association studies (e.g., Skol et al. [10]) in which only a portion of samples will be genotyped in the first stage to select markers to be genotyped using the remaining samples.

Table 1 Ranks of candidate genes among 9187 SNPs across 22 chromosomes based on five ranking methods, sorted by chromosome and location

An independent simulation study

To further study the properties of the robust ranking procedures, we conducted an independent simulation study. We simulated a case-control genome-wide association study of 100,000 SNPs with 500 cases and 500 controls. For illustration, we simply assumed that all SNPs were in linkage equilibrium, among which 9 SNPs were associated with a disease (3 SNPs had recessive, additive, and dominant modes of inheritance, respectively). The MAFs for the recessive, additive, and dominant SNPs were set at 0.3. MAFs for other null SNPs were generated from a uniform distribution (0, 1). The GRRs for each genetic model were specified. We repeated simulations of 100 K SNPs ten times and the average ranks for the 9 candidate SNPs were obtained and reported in Table 2. As in Table 1, min p and MERT are more robust than a single trend test (Z₀, Z_1/2, or Z₁) for genome-wide scans. For example, for SNPs 3, 6, and 9 (having the greatest GRRs for each genetic model), the ranks of min p and MERT across three genetic models are all on the list of top 100 most significant SNPs, but are not if any single trend test is used.

Table 2 Average ranks of nine SNPs with true association in ten replicates in a genome-wide association study with 100 K SNPs

Conclusion

In this article, we studied the robust properties of ranks of true associations in genome-wide scans. In some situations, ranking markers by a single trend test may not be robust, in particular, when the true genetic model is unknown. Using robust methods, such as min p and MERT, to rank markers may lead to higher power when the ranks by three CATTs are quite different. The results showed that they are particularly useful in ensuring that recessive effects are not missed. While min p and MERT improve the univariate approach to the first stage of gene discovery, simulated data shows that some SNPs are not found via these univariate methods.