Robust Association Tests for the Replication of Genome-Wide Association Studies

Abstract

In genome-wide association study (GWAS), robust genetic association tests such as maximum of three CATTs (MAX3), each corresponding to recessive, additive, and dominant genetic models, the minimum p value of Pearson's Chi-square test with 2 degrees of freedom, and CATT based on additive genetic model (MIN2), genetic model selection (GMS), and genetic model exclusion (GME) methods have been shown to provide better power performance under wide range of underlying genetic models. In this paper, we demonstrate how these robust tests can be applied to the replication study of GWAS and how the overall statistical significance can be evaluated using the combined test formed by p values of the discovery and replication studies.

Jungnam Joo 1 and Ju-Hyun Park 2 and Bora Lee 1 and Boram Park 1 and Sohee Kim 1 and Kyong-Ah Yoon 3 and Jin Soo Lee 3 and Nancy L. Geller 4

Academic Editor:Taesung Park

1, Biometric Research Branch, Research Institute and Hospital, National Cancer Center, Gyeonggi-do, Goyang-si 410-769, Republic of Korea
2, Department of Statistics, Dongguk University, Seoul 100-715, Republic of Korea
3, Lung Cancer Branch, Research Institute and Hospital, National Cancer Center, Gyeonggi-do, Goyang-si 410-769, Republic of Korea
4, Office of Biostatistics Research, National Heart, Lung and Blood Institute, Bethesda, MD 20892-7938, USA

Received 14 November 2014; Revised 14 February 2015; Accepted 14 February 2015; 4 August 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

With the advance of biotechnology and substantial reduction of genotyping costs, a genome-wide association study (GWAS) using hundred thousand markers in several thousand individuals is now increasingly utilized and has been successful in detecting genetic associations across the entire genome with complex human traits [1-6]. Among many challenges this application holds; development of more efficient and robust statistical methodologies with higher power to detect an association with a single marker has been one of the most important statistical issues, given that effects of individual markers are usually characterized as being small to moderate. One attempt to overcome this challenge is focused on developing efficient tests that are robust against underlying genetic model misspecification.

Two most frequently used association tests are the allele-based test (ABT) and the genotype-based test (GBT). ABT compares the allele frequencies between cases and controls, while GBT compares the genotype distributions of cases and controls. The Cochran-Armitage trend test (CATT) [7, 8] is a popular GBT which takes into account the underlying genetic model. It is well known, however, that the ABT may inflate type I error when Hardy-Weinberg equilibrium (HWE) does not hold in the samples [9]. Even under HWE, when the genetic model is recessive or dominant, the ABT may suffer from serious power loss. On the other hand, the CATT does not depend on HWE, but to apply the CATT the choice of scores optimal for the underlying genetic model needs to be specified. For complex diseases, the genetic model is usually unknown and robust tests such as the maximum of three CATTs (MAX3) [10] and the maximum efficiency robust test (MERT) [11, 12] are preferable. Alternatively, Zheng and Ng [13] and Joo et al. [14] proposed a two-phase analysis based on the genetic model selection (GMS) and genetic model exclusion (GME). Moreover, an alternative approach was proposed by the Wellcome trust case-control consortium (WTCCC) [5] which used a minimum [figure omitted; refer to PDF] value of Pearson's Chi-square test and additive CATT, and the asymptotic properties of this approach were studied in detail by Joo et al. [15]. These methods provide better or comparable power performance than some of the robust tests such as MAX3.

In this paper, we illustrate how these robust tests can be applied to a replication study of GWAS and how overall statistical significance can be evaluated using the combined test formed by [figure omitted; refer to PDF] values of the discovery and replication studies. The importance of replication or validation in GWAS has been well recognized [16, 17], and joint analysis in a two-stage design of GWAS has been proved to be more powerful than replication-based analysis and has been widely conducted in GWAS with a variety of phenotypes of interest [18, 19].

The paper is organized as follows. We first describe the data structures and notation and review existing robust association tests for a single data set. Then we describe how to obtain the [figure omitted; refer to PDF] value for the replication data set, given the significant result of the discovery stage, using robust tests. In the next section, a combined test of the [figure omitted; refer to PDF] values of the discovery and replication data sets is proposed, together with the way to evaluate the statistical significance for the combined test. Simulation studies are conducted to compare the type I error rates and powers of various analytical strategies. For illustration purposes, the summarized methods are applied to a non-small-cell lung cancer data set and at the end there is a discussion.

2. Methods

2.1. Data and Notation

For a marker with two alleles [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , let the frequencies of [figure omitted; refer to PDF] in cases and controls be [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Denote three genotypes by [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . In case-control association studies, [figure omitted; refer to PDF] cases and [figure omitted; refer to PDF] controls are independently sampled from each population. The observed genotype counts for [figure omitted; refer to PDF] are [figure omitted; refer to PDF] in the cases and [figure omitted; refer to PDF] in the controls. Disease prevalence is denoted by [figure omitted; refer to PDF] and penetrance by [figure omitted; refer to PDF] for [figure omitted; refer to PDF] . Two genotype relative risks (GRRs) are denoted by [figure omitted; refer to PDF] and [figure omitted; refer to PDF] using [figure omitted; refer to PDF] as baseline penetrance. Under the null hypothesis of no association [figure omitted; refer to PDF] or alternatively [figure omitted; refer to PDF] . Genetic model is recessive (REC), additive (ADD), multiplicative (MUL), and dominant (DOM) when [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , respectively.

2.2. Review of Association Tests for a Single Data Set

The association in case-control studies can be tested using various methods which have been extensively studied. The general association between the disease status and the SNP can be tested using Pearson's Chi-square test which has an asymptotic Chi-square distribution with 2 degrees of freedom under [figure omitted; refer to PDF] . The test is given by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] for [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Under Hardy-Weinberg equilibrium (HWE), an allele-based test (ABT) and CATT with scores [figure omitted; refer to PDF] for [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] , are given by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] [9]. The optimal choices of [figure omitted; refer to PDF] for the recessive (REC), additive/multiplicative (ADD/MUL), and dominant (DOM) models are [figure omitted; refer to PDF] and 1, respectively [9, 20]. Both [figure omitted; refer to PDF] and [figure omitted; refer to PDF] asymptotically follow a standard normal distribution under [figure omitted; refer to PDF] . [figure omitted; refer to PDF] can be used even when HWE does not hold. However, without the HWE assumption, [figure omitted; refer to PDF] does not follow a standard normal distribution due to the correlation between two alleles.

A robust test, MAX3 proposed by Friedlin et al. [10], can be obtained by taking the maximum of three CATTs under the three genetic models as [figure omitted; refer to PDF] . Parametric bootstrap or permutation methods can be used to find the [figure omitted; refer to PDF] value of [figure omitted; refer to PDF] [4].

Let the [figure omitted; refer to PDF] values of Pearson's Chi-square test and CATT under the additive genetic model [figure omitted; refer to PDF] be [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , respectively. WTCCC [5] proposed an alternative robust test [figure omitted; refer to PDF] . Joo et al. [15] derived the asymptotic null distribution of [figure omitted; refer to PDF] and using their result the [figure omitted; refer to PDF] value of [figure omitted; refer to PDF] can be obtained as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the cumulative distributions of Chi-square distributions with 1 and 2 degrees of freedom.

On the other hand, Song and Elston [21] considered a Hardy-Weinberg disequilibrium trend test (HWDTT) given by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the estimates of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Here, [figure omitted; refer to PDF] denotes the Hardy-Weinberg disequilibrium (HWD) coefficient defined by [figure omitted; refer to PDF] and [figure omitted; refer to PDF] and [figure omitted; refer to PDF] denote the HWD coefficient in cases and controls, respectively.

Zheng and Ng [13] used the information contained in the signs of [figure omitted; refer to PDF] to determine the genetic models in their two-phase method. Their two-phase statistic [figure omitted; refer to PDF] is given by [figure omitted; refer to PDF] if [figure omitted; refer to PDF] , [figure omitted; refer to PDF] if [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] otherwise, where [figure omitted; refer to PDF] for [figure omitted; refer to PDF] . The asymptotic correlations between [figure omitted; refer to PDF] and three CATTs under HWE were derived and the significance level was adjusted accordingly to control the desired type I error. Based on the observation that this method assumes [figure omitted; refer to PDF] is the risk allele, Joo et al. [14] studied the behavior of [figure omitted; refer to PDF] when either one of the alleles can be a risk allele. They chose the risk allele based on the sign of [figure omitted; refer to PDF] ; that is, if [figure omitted; refer to PDF] , [figure omitted; refer to PDF] is the risk allele, and [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are chosen for REC, ADD, and DOM models, respectively. If [figure omitted; refer to PDF] , the respective test statistics are chosen to be [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . They incorporate this property in defining the test statistic for genetic model selection ( [figure omitted; refer to PDF] ) and calculating the [figure omitted; refer to PDF] value. Let [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Then, the [figure omitted; refer to PDF] value of this method can be obtained by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] in (5) and [figure omitted; refer to PDF] ( [figure omitted; refer to PDF] ) are replaced by their consistent estimates. Here, [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Moreover, [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the observed values of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , respectively.

While studying the properties of GMS, Joo et al. [14] noticed that the probability of selecting the true recessive or dominant models using [figure omitted; refer to PDF] is very low especially for low to moderate GRRs, but the unlikely genetic model can be successfully excluded. This led to genetic model exclusion method [figure omitted; refer to PDF] which is the same as the [figure omitted; refer to PDF] described above except [figure omitted; refer to PDF] for [figure omitted; refer to PDF] is replaced by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] . And the [figure omitted; refer to PDF] value of GME can be obtained as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] for [figure omitted; refer to PDF] .

2.3. [figure omitted; refer to PDF] Value of Replication Data Using the Robust Method

In the discovery stage, the [figure omitted; refer to PDF] value of robust association tests, including [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , can be obtained as described in Section 2.2. For the [figure omitted; refer to PDF] value of replication data using the robust method, we use the same analytic method that was used for discovery and the risk allele identified by it [16]. This means that when the best test statistic or genetic model is selected in the discovery stage, the replication stage will adopt the discovery stage selection and the direction of association.

Suppose that, for simplicity of notation, our interest is in GWAS with two stages, one for discovery and the other for replication, although the methodology described below can be extended to multistages for replication. Let [figure omitted; refer to PDF] for [figure omitted; refer to PDF] be the CATT optimal for recessive, additive, and dominant models and let [figure omitted; refer to PDF] be corresponding [figure omitted; refer to PDF] value for [figure omitted; refer to PDF] th stage ( [figure omitted; refer to PDF] for discovery and [figure omitted; refer to PDF] for replication stages). Also, denote [figure omitted; refer to PDF] for [figure omitted; refer to PDF] . Then, for CATT with a preselected genetic model, [figure omitted; refer to PDF] using a one-sided [figure omitted; refer to PDF] value given the direction of association from the discovery stage, and [figure omitted; refer to PDF] . Moreover, denote the test statistics and [figure omitted; refer to PDF] values using Pearson's Chi-square test from the [figure omitted; refer to PDF] th stage as [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Further, let HWDTT from the [figure omitted; refer to PDF] th stage be [figure omitted; refer to PDF] . Then, the second stage [figure omitted; refer to PDF] values, using [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , denoted as [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , can be obtained as follows: [figure omitted; refer to PDF]

It is important to note that even though the direction of the test statistics and the selected genetic models are used to obtain the second stage [figure omitted; refer to PDF] values, the [figure omitted; refer to PDF] values from the two stages are independent under the null hypothesis. This is because, under the null hypothesis, the probability of [figure omitted; refer to PDF] being positive or negative is simply [figure omitted; refer to PDF] , and the probability of the selection of a certain genetic model is also a constant ( [figure omitted; refer to PDF] for the recessive and dominant models and [figure omitted; refer to PDF] for the additive model).

2.4. Combined Test Using [figure omitted; refer to PDF] Values and Its Statistical Significance

For a given robust test, we can consider the joint analysis by combining [figure omitted; refer to PDF] values from the discovery and replication stages of GWAS. We consider using [figure omitted; refer to PDF] values rather than the test statistics because test statistics can have complex forms and obtaining the distribution of the joint test can be difficult. On the other hand, calculating a [figure omitted; refer to PDF] value for each data set might be relatively simple, and the distribution of [figure omitted; refer to PDF] values under the null hypothesis of no association is easy to handle.

There are several methods for combining test statistics from two stages [22], and two most commonly used forms are based on Fisher's combination and a linear combination after inverse normal transformation [23]. Fisher's combination (FC) directly sums [figure omitted; refer to PDF] values after [figure omitted; refer to PDF] transformation; that is, [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] is [figure omitted; refer to PDF] value from [figure omitted; refer to PDF] for discovery and [figure omitted; refer to PDF] for replication stages using a given robust test. A specification of [figure omitted; refer to PDF] gives the same weight for discovery and replication stages, and one can consider [figure omitted; refer to PDF] and [figure omitted; refer to PDF] where [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are sample sizes of the discovery and replication data sets. A linear combination of two [figure omitted; refer to PDF] values after taking the inverse of the standard normal cumulative distribution is given by [figure omitted; refer to PDF] with a natural choice of [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . Let the significance level of the discovery stage be [figure omitted; refer to PDF] , which means that markers with [figure omitted; refer to PDF] are selected and replicated in the replication stage. The [figure omitted; refer to PDF] value of combined test can then be obtained as [figure omitted; refer to PDF] where the observed value of [figure omitted; refer to PDF] is [figure omitted; refer to PDF] . The [figure omitted; refer to PDF] are calculated as [figure omitted; refer to PDF] for equal weights where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] for unequal weights where [figure omitted; refer to PDF] . Detailed derivations are described in the Appendix. Equivalently, for an overall type I error threshold for a single marker of [figure omitted; refer to PDF] , one may obtain the threshold [figure omitted; refer to PDF] of [figure omitted; refer to PDF] that satisfies [figure omitted; refer to PDF] . Similarly, for the [figure omitted; refer to PDF] , the [figure omitted; refer to PDF] value is calculated as [figure omitted; refer to PDF] = [figure omitted; refer to PDF] for [figure omitted; refer to PDF] where the observed value of [figure omitted; refer to PDF] .

3. Simulation Results

3.1. Type I Error

Table 1 provides the type I errors under different scenarios. A disease prevalence of 10% is assumed, and a total of 1500 cases and 1500 controls were divided into two stages. The proportions of samples in the first stage ( [figure omitted; refer to PDF] ) of 0.5 and 0.6 were considered for the minor allele frequency (MAF) of 0.3 and 0.4. We considered [figure omitted; refer to PDF] markers to control the genome-wide false positive rate at [figure omitted; refer to PDF] with the Bonferroni correction. We did not consider a larger [figure omitted; refer to PDF] such as 300,000 or 500,000 because this will require more than 10 million simulations to show a stable estimate of the type I error rate. With [figure omitted; refer to PDF] , we performed 20,000 simulations which result in less than 10% of a coefficient of variation for a significance level [figure omitted; refer to PDF] for each marker [24]. The test statistics considered are [figure omitted; refer to PDF] , Pearson's Chi-square test, MIN2, MAX3, GMS, and GME. For the second stage analysis, we considered a replication-based analysis, [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] as proposed above. The results are based on the situation under HWE (HWE coefficient [figure omitted; refer to PDF] ). As expected, all tests control the type I error reasonablly well, and similar results were obtained when a slight deviation from HWE is present with [figure omitted; refer to PDF] (results not shown).

Table 1: Type I error rates of three approaches-replication-based (REP) test, Fisher's combination ( [figure omitted; refer to PDF] ), and linear combination of test ( [figure omitted; refer to PDF] )-based on the CATT with an additive model ( [figure omitted; refer to PDF] ), [figure omitted; refer to PDF] , MAX3, MIN2, GMS, and GME. The disease prevalence [figure omitted; refer to PDF] , [figure omitted; refer to PDF] markers, [figure omitted; refer to PDF] cases, and [figure omitted; refer to PDF] controls are considered based on 20,000 simulations.

MAF	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]		[figure omitted; refer to PDF]
[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	MAX3	MIN2	GMS	GME
0.3	0.5	0.05	REP	0.00530	0.00455	0.00505	0.00485	0.0050	0.00490
[figure omitted; refer to PDF]	0.00500	0.00535	0.00495	0.00510	0.00515	0.00460
[figure omitted; refer to PDF]	0.00535	0.00525	0.00485	0.00510	0.0050	0.00485
0.1	REP	0.00510	0.00560	0.00565	0.00485	0.00525	0.00545
[figure omitted; refer to PDF]	0.00565	0.00535	0.00565	0.00545	0.00565	0.00540
[figure omitted; refer to PDF]	0.00520	0.00565	0.00525	0.00520	0.00530	0.00525
0.3	0.6	0.05	REP	0.00510	0.00485	0.00480	0.00515	0.00480	0.00500
[figure omitted; refer to PDF]	0.00445	0.00455	0.00450	0.00455	0.00450	0.00460
LC	0.00500	0.00515	0.00495	0.00520	0.00475	0.00480
0.1	REP	0.00500	0.00485	0.00485	0.00535	0.00530	0.00505
[figure omitted; refer to PDF]	0.00465	0.00490	0.00485	0.00500	0.00455	0.00460
[figure omitted; refer to PDF]	0.00480	0.00470	0.00515	0.00490	0.00485	0.00475
0.4	0.5	0.05	REP	0.00590	0.00505	0.00530	0.00565	0.00505	0.00510
[figure omitted; refer to PDF]	0.00575	0.00430	0.00460	0.00535	0.00460	0.00500
[figure omitted; refer to PDF]	0.00600	0.00445	0.00500	0.00540	0.00490	0.00490
0.1	REP	0.00525	0.00470	0.00535	0.00450	0.00480	0.00515
[figure omitted; refer to PDF]	0.00515	0.00510	0.00495	0.00475	0.00540	0.00500
[figure omitted; refer to PDF]	0.00530	0.00500	0.00485	0.00475	0.00495	0.00510
0.4	0.6	0.05	REP	0.00475	0.00585	0.00480	0.00500	0.00515	0.00495
[figure omitted; refer to PDF]	0.00460	0.00470	0.00420	0.00490	0.00455	0.00440
[figure omitted; refer to PDF]	0.00525	0.00550	0.00520	0.00580	0.00510	0.00510
0.1	REP	0.00550	0.00490	0.00515	0.00535	0.00555	0.00540
[figure omitted; refer to PDF]	0.00520	0.00370	0.00495	0.00450	0.00515	0.00510
[figure omitted; refer to PDF]	0.00565	0.00485	0.00570	0.00530	0.00610	0.00580

3.2. Empirical Power

We examined the empirical powers of different tests considered above. In Figure 1, we considered [figure omitted; refer to PDF] markers, a disease prevalence of 10%, the same genotype relative risk for two stages ( [figure omitted; refer to PDF] and [figure omitted; refer to PDF] ), and 1,000 cases and 1,000 controls. 2,000 simulations were performed under HWE ( [figure omitted; refer to PDF] ) to control the genome-wide false positive rate at [figure omitted; refer to PDF] . The recessive, additive, and dominant models were assumed for the first, second, and third rows. Both joint analyses showed better power performances compared to the replication-based analysis (up to 15.9% in scenarios considered in Figure 1), and LC and FC have comparable powers with less than 2% difference. The power gain of using the joint analysis is not as much as that observed in Skol et al. [18]. However, as reported by Skol et al. [18], when the between-stage heterogeneity exists and the risk allele has a larger effect in the first stage than that in the second stage, much improved power is observed by using the joint test. Figure 2 shows results under this scenario with [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , and the observed increase in power using the joint test is as high as 33.9%. Again, the difference between LC and FC is minor with less than 3% difference. As for comparison between different robust methods, MAX3, GMS, and GME perform well under the recessive model, while [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and MIN2 are less powerful. Under the additive model, [figure omitted; refer to PDF] is most powerful, as expected, and [figure omitted; refer to PDF] is least powerful. Other robust methods perform well with a slight decrease in power compared to [figure omitted; refer to PDF] . Under the dominant model, MAX3, GMS, and GME perform the best even though all tests show good power performances, and the difference is minor. Similar patterns were observed when a slight deviation from the HWE is present (results not shown).

Figure 1: Empirical powers based on 2,000 simulations for [figure omitted; refer to PDF] markers, genotype relative risks of both stages = 1.4, and disease prevalence [figure omitted; refer to PDF] under the recessive, additive, and dominant models. 1,000 cases and 1,000 controls are considered to control [figure omitted; refer to PDF] . The first stage type I error rate for discovery is [figure omitted; refer to PDF] . Six test statistics, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , MAX3, MIN2, GMS, and GME, are considered. The first, second, and third columns depict powers using the replication-based test, [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , respectively.

[figure omitted; refer to PDF]

Figure 2: Empirical powers based on 2,000 simulations for [figure omitted; refer to PDF] markers; genotype relative risks of two stages are different ( [figure omitted; refer to PDF] , [figure omitted; refer to PDF] ); disease prevalence [figure omitted; refer to PDF] under the recessive, additive, and dominant models. 1,000 cases and 1,000 controls are considered to control [figure omitted; refer to PDF] . The first stage type I error rate for discovery is [figure omitted; refer to PDF] . Six test statistics, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , MAX3, MIN2, GMS, and GME, are considered. The first, second, and third columns depict powers using the replication-based test, [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] , respectively.

[figure omitted; refer to PDF]

4. Real Data Application

The GWAS on non-small-cell lung cancer (NSCLC) by Yoon et al. [25] studied 621 NSCLC patients and 1541 control subjects in the discovery stage. After stringent quality control steps, a total of 246,758 SNPs were tested for the association with NSCLC based on [figure omitted; refer to PDF] . In the replication stage, 168 SNPs with [figure omitted; refer to PDF] value less than [figure omitted; refer to PDF] in the first stage based on [figure omitted; refer to PDF] were tested using 804 patients and 1470 control samples. We identified additional 234 SNPs using MIN2 in the first stage which could be studied in the replication stage if MIN2 was used instead of [figure omitted; refer to PDF] since MIN2 produces stronger evidence for the additional SNPs than [figure omitted; refer to PDF] does. The Manhattan plots of using MIN2 and [figure omitted; refer to PDF] are presented in Figure 3. One example is [figure omitted; refer to PDF] located in chromosome 2, which had a [figure omitted; refer to PDF] value of [figure omitted; refer to PDF] which reached significance level at Bonferroni correction in discovery samples alone, whereas [figure omitted; refer to PDF] yielded a [figure omitted; refer to PDF] value greater than [figure omitted; refer to PDF] . Even though there is possibility of false positive findings, these SNPs could have been selected for replication if robust methods were used.

Figure 3: Manhattan plots of 246,758 SNPs from Yoon et al. [25] based on MIN2 (a) and [figure omitted; refer to PDF] (b). The [figure omitted; refer to PDF] axis is chromosomal location and the [figure omitted; refer to PDF] axis is the significance ( [figure omitted; refer to PDF] ) of association. The horizontal line corresponds to the significance level [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Since we do not have replication data for these additional SNPs selected using MIN2 because the first stage selection was based on [figure omitted; refer to PDF] in Yoon et al. [25], just for illustration purpose of the proposed methods, we present the results of three SNPs including [figure omitted; refer to PDF] that was reported by Yoon et al. [25]. When the significance level in the discovery stage is set at [figure omitted; refer to PDF] so that all these exemplary SNPs can be selected in the discovery stage; the [figure omitted; refer to PDF] value of combined test based on four robust methods ( [figure omitted; refer to PDF] , MIN2, GMS, and GME) as well as [figure omitted; refer to PDF] and Pearson's Chi-square test is presented in Table 2. Fisher's combination was used for the joint test in the second stage. Only [figure omitted; refer to PDF] was found to be significant with Bonferroni correction ( [figure omitted; refer to PDF] value [figure omitted; refer to PDF] ) by all except [figure omitted; refer to PDF] method.

Table 2: For selected exemplary three SNPs for testing association with NSCLC, [figure omitted; refer to PDF] value of combined test using additive CATT ( [figure omitted; refer to PDF] ), Pearson's Chi-square test ( [figure omitted; refer to PDF] ), MAX3, MIN2, [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] .

SNP	[figure omitted; refer to PDF] value of [figure omitted; refer to PDF]	[figure omitted; refer to PDF] value of [figure omitted; refer to PDF]
	Discovery	Replication	Combined test	Discovery	Replication	Combined test
rs2131877	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
rs905551	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
rs1695109	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
SNP	[figure omitted; refer to PDF] value of MAX3	[figure omitted; refer to PDF] value of MIN2
Discovery	Replication	Combined test	Discovery	Replication	Combined test
rs2131877	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
rs905551	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
rs1695109	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
SNP	[figure omitted; refer to PDF] value of [figure omitted; refer to PDF]	[figure omitted; refer to PDF] value of [figure omitted; refer to PDF]
Discovery	Replication	Combined test	Discovery	Replication	Combined test
rs2131877	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
rs905551	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]
rs1695109	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]	[figure omitted; refer to PDF]

5. Discussion

In genetic association studies, efficiency robust tests whose performance does not depend on the underlying genetic model have been extensively studied, and their power benefit over a wide range of genetic models has been well recognized. In this paper, we described how the idea of these robust association tests can be applied to the replication studies and further how overall statistical significance can be evaluated using the combined test formed by [figure omitted; refer to PDF] values of the discovery and replication studies.

When the robust tests are used, the test statistic of each stage can have a complex form and thus dealing with the distribution of the joint test can be difficult, whereas calculating the [figure omitted; refer to PDF] value of each stage might be relatively simple. Because the asymptotic distribution of the [figure omitted; refer to PDF] value under the null hypothesis of no association is easy to handle, the combined test using [figure omitted; refer to PDF] values rather than the test statistics themselves can provide computational convenience.

There are several methods for combining test statistics from two stages and Won et al. [22] compared the performances of various choices. Two most commonly used forms are based on Fisher's combination and the linear combination after the inverse normal transformation [23], and we presented the test statistics and [figure omitted; refer to PDF] values of these two methods. In our limited experience, the linear combination and Fisher's combination are fairly comparable. Fisher's combination seems to perform slightly better than the linear combination when there exists some heterogeneity between stages in terms of the genotype relative risk, while the linear combination seems to perform slightly better in most of other situations. However, the difference is extremely minor. Further research is required for the thorough comparison of various methods of combining [figure omitted; refer to PDF] values in the application of efficiency robust tests to the replication of genetic association studies.

In a genetic study where the purpose of considering a replication stage is to validate or replicate the genetic findings from the discovery stage, which is the case considered in this paper, the analysis in the replication stage utilized the test statistic or genetic model that is selected as being the best in the discovery stage and also the direction of the risk allele, following guidelines for exact replication in genetic association studies. If the purpose is to simply combine the evidence from different data sources such as in meta-analysis, other strategies may be devised. Further research, again, is required to provide fully detailed properties of such methods.

Power gain of a joint analysis over the conventional replication-based analysis was thoroughly studied by Skol et al. [18, 19]. In our simulation, the amount of power increase using a joint test compared to the replication-based analysis was much minor than what was observed by Skol et al. [18, 19]. The exact reason is not known, but we suspect this might be due to the power advantages of robust methods and also due to the fact that the optimal choice from the first stage is used when calculating the second stage [figure omitted; refer to PDF] values. However, even though it was minor in some situations, the joint anlysis presented better power performance than the replication-based analysis in our study. This type of joint analysis raised concerns about the exact meaning of replication [17]. However, McCarthy et al. [26] mentioned that joint analyses "blur the boundaries of where exactly replication starts, but whichever analytical approach is taken, confirmation in many independent samples is important and it is the overall strength of the evidence of association that matters." Purpose of the current study was to present how the overall strength of the evidence of association can be evaluated when robust tests are used in GWAS replication studies.

We illustrated how the proposed methods can be applied in the real data that studied the association of SNPs with non-small-cell lung cancer (NSCLC) in discovery and replication stages. In the original study reported by Yoon et al. [25], SNPs were selected in the discovery data set not based on the robust tests but based on additive CATT. Therefore, we found that some SNPs could have been selected by one of the robust methods but they were not included in the replication data set. For these SNPs, we were not able to perform the joint analysis that we propose, and it was not possible to examine whether there are other SNPs that could have been found to be associated with NSCLC by proposed methods in the replication study. For this reason, we merely presented how many additional SNPs could have been further followed in the replication stage when robust methods were used. In many GWASs, it is a common practice to report the summary test statistics and [figure omitted; refer to PDF] values of the SNPs under a specific genetic model, usually an additive model, which were further genotyped in the replication stage and were finally defined to be significantly associated with a phenotype of interest. As emphasized in this paper, one may have a better chance of finding many missing SNPs by applying more powerful and robust methods that consider different genetic models simultaneously. Therefore, we urge the community to share test results under not only an additive model but also other genetic models, although they were not significant at a stringent significance level, so that future research may have enriched data resources, to which robust tests can be applied in association studies.

Acknowledgments

The authors are indebted to late Dr. Gang Zheng for his inspiration and support on their work. This work was supported by grant of the National Cancer Center (no. NCC-1210060).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] R. J. Klein, C. Zeiss, E. Y. Chew, J.-Y. Tsai, R. S. Sackler, C. Haynes, A. K. Henning, J. P. SanGiovanni, S. M. Mane, S. T. Mayne, M. B. Bracken, F. L. Ferris, J. Ott, C. Barnstable, J. Hoh, "Complement factor H polymorphism in age-related macular degeneration," Science , vol. 308, no. 5720, pp. 385-389, 2005.

[2] R. H. Duerr, K. D. Taylor, S. R. Brant, J. D. Rioux, M. S. Silverberg, M. J. Daly, A. H. Steinhart, C. Abraham, M. Regueiro, A. Griffiths, T. Dassopoulos, A. Bitton, H. Yang, S. Targan, L. W. Datta, E. O. Kistner, L. P. Schumm, A. T. Lee, P. K. Gregersen, M. M. Barmada, J. I. Rotter, D. L. Nicolae, J. H. Cho, "A genome-wide association study identifies IL 23R as an inflammatory bowel disease gene," Science , vol. 314, no. 5804, pp. 1461-1463, 2006.

[3] A. Herbert, N. P. Gerry, M. B. McQueen, I. M. Heid, A. Pfeufer, T. Illig, H.-E. Wichmann, T. Meitinger, D. Hunter, F. B. Hu, G. Colditz, A. Hinney, J. Hebebrand, K. Koberwitz, X. Zhu, R. Cooper, K. Ardlie, H. Lyon, J. N. Hirschhorn, N. M. Laird, M. E. Lenburg, C. Lange, M. F. Christman, "A common genetic variant is associated with adult and childhood obesity," Science , vol. 312, no. 5771, pp. 279-283, 2006.

[4] R. Sladek, G. Rocheleau, J. Rung, C. Dina, L. Shen, D. Serre, P. Boutin, D. Vincent, A. Belisle, S. Hadjadj, B. Balkau, B. Heude, G. Charpentier, T. J. Hudson, A. Montpetit, A. V. Pshezhetsky, M. Prentki, B. I. Posner, D. J. Balding, D. Meyre, C. Polychronakos, P. Froguel, "A genome-wide association study identifies novel risk loci for type 2 diabetes," Nature , vol. 445, no. 7130, pp. 881-885, 2007.

[5] P. R. Burton, D. G. Clayton, L. R. Cardon, N. Craddock, P. Deloukas, A. Duncanson, D. P. Kwiatkowski, M. I. McCarthy, W. H. Ouwehand, N. J. Samani, J. A. Todd, P. Donnelly, J. C. Barrett, D. Davison, D. Easton, D. Evans, H.-T. Leung, J. L. Marchini, A. P. Morris, C. C. A. Spencer, M. D. Tobin, "Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls," Nature , vol. 447, no. 7145, pp. 661-678, 2007.

[6] S. H. Kwak, S. H. Kim, Y. M. Cho, M. J. Go, Y. S. Cho, S. H. Choi, M. K. Moon, H. S. Jung, H. D. Shin, H. M. Kang, N. H. Cho, I. K. Lee, S. Y. Kim, B.-G. Han, H. C. Jang, K. S. Park, "A genome-wide association study of gestational diabetes mellitus in Korean women," Diabetes , vol. 61, no. 2, pp. 531-541, 2012.

[7] W. G. Cochran, "Some methods for strengthening the common ^χ2 tests," Biometrics , vol. 10, no. 4, pp. 417-451, 1954.

[8] P. Armitage, "Tests for linear trends in proportions and frequencies," Biometrics , vol. 11, no. 3, pp. 375-386, 1955.

[9] P. D. Sasieni, "From genotypes to genes: doubling the sample size," Biometrics , vol. 53, no. 4, pp. 1253-1261, 1997.

[10] B. Freidlin, G. Zheng, Z. Li, J. L. Gastwirth, "Trend tests for case-control studies of genetic markers: power, sample size and robustness," Human Heredity , vol. 53, no. 3, pp. 146-152, 2002.

[11] J. L. Gastwirth, "On robust procedures," Journal of the American Statistical Association , vol. 61, no. 316, pp. 929-948, 1966.

[12] J. L. Gastwirth, "The use of maximin efficiency robust tests in combining contingency tables and survival analysis," Journal of the American Statistical Association , vol. 80, no. 390, pp. 380-384, 1985.

[13] G. Zheng, H. K. T. Ng, "Genetic model selection in two-phase analysis for case-control association studies," Biostatistics , vol. 9, no. 3, pp. 391-399, 2008.

[14] J. Joo, M. Kwak, G. Zheng, "Improving power for testing genetic association in case-control studies by reducing the alternative space," Biometrics , vol. 66, no. 1, pp. 266-276, 2010.

[15] J. Joo, M. Kwak, K. Ahn, G. Zheng, "A robust genome-wide scan statistic of the Wellcome Trust Case-Control Consortium," Biometrics , vol. 65, no. 4, pp. 1115-1122, 2009.

[16] P. Kraft, E. Zeggini, J. P. A. Ioannidis, "Replication in genome-wide association studies," Statistical Science , vol. 24, no. 4, pp. 561-573, 2009.

[17] D. C. Thomas, G. Casey, D. V. Conti, R. W. Haile, J. P. Lewinger, D. O. Stram, "Methodological issues in multistage genome-wide association studies," Statistical Science , vol. 24, no. 4, pp. 414-429, 2009.

[18] A. D. Skol, L. J. Scott, G. R. Abecasis, M. Boehnke, "Joint analysis is more efficient than replication-based analysis for two-stage genome-wide association studies," Nature Genetics , vol. 38, no. 2, pp. 209-213, 2006.

[19] A. D. Skol, L. J. Scott, G. R. Abecasis, M. Boehnke, "Optimal designs for two-stage genome-wide association studies," Genetic Epidemiology , vol. 31, no. 7, pp. 776-788, 2007.

[20] G. Zheng, B. Freidlin, Z. Li, J. L. Gastwirth, "Choice of scores in trend tests for case-control studies of candidate-gene associations," Biometrical Journal , vol. 45, no. 3, pp. 335-348, 2003.

[21] K. Song, R. C. Elston, "A powerful method of combining measures of association and Hardy-Weinberg Disequilibrium for fine-mapping in case-control studies," Statistics in Medicine , vol. 25, no. 1, pp. 105-126, 2006.

[22] S. Won, N. Morris, Q. Lu, R. C. Elston, "Choosing an optimal method to combine P -values," Statistics in Medicine , vol. 28, no. 11, pp. 1537-1553, 2009.

[23] F. Begum, D. Ghosh, G. C. Tseng, E. Feingold, "Comprehensive literature review and statistical considerations for GWAS meta-analysis," Nucleic Acids Research , vol. 40, no. 9, pp. 3777-3784, 2012.

[24] B. Efron, R. J. Tibshirani An Introduction to the Bootstrap , Chapman & Hall/CRC Press, 1993.

[25] K. A. Yoon, J. H. Park, J. Han, S. Park, G. K. Lee, J.-Y. Han, J. I. Zo, J. Kim, J. E. Lee, A. Takahashi, M. Kubo, Y. Nakamura, J. S. Lee, "A genome-wide association study reveals susceptibility variants for non-small cell lung cancer in the Korean population," Human Molecular Genetics , vol. 19, no. 24, pp. 4948-4954, 2010.

[26] M. I. McCarthy, G. R. Abecasis, L. R. Cardon, D. B. Goldstein, J. Little, J. P. A. Ioannidis, J. N. Hirschhorn, "Genome-wide association studies for complex traits: consensus, uncertainty and challenges," Nature Reviews Genetics , vol. 9, no. 5, pp. 356-369, 2008.

[27] I. J. Good, "On the weighted combination of significance tests," Journal of the Royal Statistical Society, Series B , vol. 17, pp. 264-265, 1955.

Appendix

Appendix

[figure omitted; refer to PDF] value of Fisher's Combination for Equal and Unequal Weights

Equal Weights . Assume [figure omitted; refer to PDF] . Under the null hypothesis of no association, [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are independent and each asymptotically follows a [figure omitted; refer to PDF] distribution with 2 degrees of freedom Let [figure omitted; refer to PDF] and [figure omitted; refer to PDF] be the probability and cumulative density functions of [figure omitted; refer to PDF] random variable with [figure omitted; refer to PDF] degrees of freedom Then [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . Denote the cutoff of the discovery stage based on [figure omitted; refer to PDF] as [figure omitted; refer to PDF] ; that is, [figure omitted; refer to PDF] . For observed value [figure omitted; refer to PDF] of [figure omitted; refer to PDF] , the [figure omitted; refer to PDF] value is written as [figure omitted; refer to PDF]

Unequal Weights. When different proportions of samples are used in the discovery and replication stages, it may be more appropriate to assign weights proportional to the sample sizes for each stage. For example, when only a small portion is used in the discovery stage, to prevent Fisher's combination test from being dominated by the significant result in the discovery stage, one may want to assign a small weight to the discovery stage result.

When [figure omitted; refer to PDF] is the proportion of samples used in the discovery stage, one selection for weights is [figure omitted; refer to PDF] and [figure omitted; refer to PDF] for discovery and replication stages, which simplifies to equal weights when [figure omitted; refer to PDF] . Based on these weights, we consider unequal-weighted Fisher's combination as [figure omitted; refer to PDF] [27]. Its density function is given by [figure omitted; refer to PDF] and the probability distribution function is [figure omitted; refer to PDF]

Using the previous notation, we have the following form of [figure omitted; refer to PDF] value: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] .