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

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 misspecification. 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 comfrom Genetic Analysis Workshop 15 St. Pete Beach, Florida, USA. 11–15 November 2006 Published: 18 December 2007 BMC Proceedings 2007, 1(Suppl 1):S165 <supplement> <title> <p>Genetic Analysis Workshop 15: Gene Expression Analysis and Approaches to Detecting Multiple Functional Loci</p> </title> <editor>Heather J Cordell, Mariza de Andrade, Marie-Claude Babron, Christopher W Bartlett, Joseph Beyene, Heike Bickeböller, Robert Culverhouse, Adrienne Cupples, E Warwick Daw, Josée Dupuis, Catherine T Falk, Saurabh Ghosh, Katrina A Goddard, Ellen L Goode, Elizabeth R Hauser, Lisa J Martin, Maria Martinez, Kari E North, Nancy L Saccone, Silke Schmidt, William Tapper, Duncan Thomas, David Tritchler, Veronica J Vieland, Ellen M Wijsman, Marsha A Wilcox, John S Witte, Qio g Yang, Andreas Ziegler, Laura Almasy a d Jean W MacCluer</editor> <note>Proce dings</note> <url>http://www.biomedc ntral.com/content/pdf/1753-6561-1-S1-info.pdf</url> </supplement> This article is available from: http://www.biomedcentral.com/1753-6561/1/S1/S165 © 2007 Zheng et al; licensee BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 probabil-ity 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 com-plex 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.

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 0 = dd, G 1 = Dd, and G 2 = DD are denoted as (r 0 , r 1 , r 2 ) in cases and (s 0 , s 1 , s 2 ) in controls, which follow multinomial distributions mul(r: p 0 , p 1 , p 2 ) and mul(s: q 0 , q 1 , q 2 ), 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, Without loss of generality, assume that D has high risk. Then the null hypothesis of no association can be stated as H 0 : f 0 = f 1 = f 2 = k. The alternative hypothesis is H 1 : f 0 ≤ f 1 ≤ f 2 with at least one inequality. The genotype relative risks (GRRs) are defined as λ 1 = f 1 /f 0 and λ 2 = f 2 /f 0 . The recessive, additive, and dominant models are referred to as λ 1 = 1, λ 1 = (1 + λ 2 )/2, and λ 1 = λ 2 , respectively [2][3][4].

Ranking markers with multiple statistics
When the genetic model is unknown, the three CATTs (Z 0 , Z 1/2 , Z 2 ) 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.

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 0 (Z 1/2 and Z 1 ), 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 Z n x sr rs rsn n x n x n 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 1 , and 8172 when Z 0 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.

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 0 , Z 1/2 , or Z 1 ) for genome-wide scans. For example, for SNPs 3, 6, and 9 (having the greatest GRRs for each genetic

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.