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

Authors: Jungnam Joo; Ju-Hyun Park; Bora Lee; Boram Park; Sohee Kim; Kyong-Ah Yoon; Jin Soo Lee; Nancy L. Geller
Journal: BioMed Research International
 (2015)

---

## 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 minimump 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.


---

## Body

## 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 p 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 byp 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 thep 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 p 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 allelesA and B, let the frequencies of B in cases and controls be p=P(B∣case) and q=P(B∣control). Denote three genotypes by G0=AA, G1=AB, and G2=BB. In case-control association studies, r cases and s controls are independently sampled from each population. The observed genotype counts for (G0,G1,G2) are (r0,r1,r2) in the cases and (s0,s1,s2) in the controls. Disease prevalence is denoted by k=P(disease) and penetrance by fi=P(disease∣Gi) for i=0,1,2. Two genotype relative risks (GRRs) are denoted by λ1=f1/f0 and λ2=f2/f0 using f0>0 as baseline penetrance. Under the null hypothesis of no association H0:f0=f1=f2=k or alternatively H0:λ2=λ1=1. Genetic model is recessive (REC), additive (ADD), multiplicative (MUL), and dominant (DOM) when λ1=1, λ1=(1+λ2)/2, λ1=λ21/2, and λ2=λ1, 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 underH0. The test is given by (1)Tchi2=∑j=02‍rj-njr/n2njr/n+∑j=02‍sj-njs/n2njs/n,where ni=ri+si for i=0,1,2 and n=r+s. Under Hardy-Weinberg equilibrium (HWE), an allele-based test (ABT) and CATT with scores (0,x,1) for (G0,G1,G2), where 0≤x≤1, are given by (2)ZABT=n1/22r2s0+s1-2s2r0+r12rs2n0+n1n1+2n21/2,Zx=n1/2∑i=02xisri-rsirsnn∑i=02xi2ni-∑i=02xini21/2,where (x0,x1,x2)=(0,x,1) [9]. The optimal choices of x for the recessive (REC), additive/multiplicative (ADD/MUL), and dominant (DOM) models are x=0,1/2 and 1, respectively [9, 20]. Both Zx and ZABT asymptotically follow a standard normal distribution under H0. Zx can be used even when HWE does not hold. However, without the HWE assumption, ZABT 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 MAX3=max⁡(|Z0|,|Z1/2|,|Z1|). Parametric bootstrap or permutation methods can be used to find the p value of MAX3 [4].Let thep values of Pearson’s Chi-square test and CATT under the additive genetic model Z1/2 be Pchi2 and P1/2, respectively. WTCCC [5] proposed an alternative robust test MIN2=min⁡(Pchi2,P1/2). Joo et al. [15] derived the asymptotic null distribution of MIN2 and using their result the p value of MIN2 can be obtained as (3)PMIN2=12exp⁡-12H1-1(1-MIN2)+12MIN2-12π∫H1-11-MIN2-2log⁡⁡MIN2‍e-v/2arcsin2H1-11-MIN2v-1dv,where H1 and H2 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 (4)ZH=rs/n1/2(Δ^p-Δ^q)1-n2/n-n1/(2n)n2/n+n1/(2n),where Δ^p=p^2-(p^2+p^1/2)2 and Δ^q=q^2-(q^2+q^1/2)2 are the estimates of Δp and Δq, where p^i=ri/r and q^i=si/s. Here, Δ denotes the Hardy-Weinberg disequilibrium (HWD) coefficient defined by Pr(BB)-{Pr(AB)/2+Pr(BB)}2 and Δp and Δq denote the HWD coefficient in cases and controls, respectively.Zheng and Ng [13] used the information contained in the signs of (Δp,Δq) to determine the genetic models in their two-phase method. Their two-phase statistic ZGMS is given by ZGMS=Z0 if ZH>c, Z1 if ZH<-c, and Z1/2 otherwise, where c=Φ-1(1-αH) for αH=0.05. The asymptotic correlations between ZH 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 B is the risk allele, Joo et al. [14] studied the behavior of ZGMS when either one of the alleles can be a risk allele. They chose the risk allele based on the sign of Z1/2; that is, if Z1/2>0, B is the risk allele, and Z0, Z1/2, and Z1 are chosen for REC, ADD, and DOM models, respectively. If Z1/2<0, the respective test statistics are chosen to be -Z1, -Z1/2, and -Z0. They incorporate this property in defining the test statistic for genetic model selection (ZGMS) and calculating the p value. Let Θ0(z)={z:z>c}, Θ1/2(z)={z:|z|<c}, and Θ1(z)={z:z<-c}. Then, the p value of this method can be obtained by(5)pGMS=2∑x=01‍∫Θx(zH)‍∫0∞‍∫t∞‍ϕx(zx,z1/2,zH)dzxdz1/2dzH+2∫Θ1/2(z)‍Φ(-t∧0)+ρ1/2z1-ρ1/221/2dΦ(z),where ρx=Corr(Zx,ZH) in (5) and ρx,1/2=Corr(Zx,Z1/2) (x=0,1) are replaced by their consistent estimates. Here, t=zGMS and (-t∧0)=min⁡(-t,0). Moreover, zGMS and z1/2 are the observed values of ZGMS and Z1/2, respectively.While studying the properties of GMS, Joo et al. [14] noticed that the probability of selecting the true recessive or dominant models using ZH 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 ZGME which is the same as the ZGMS described above except Zx for x=0,1/2,1 is replaced by Zx* where Zx*=(Zx+Z1/2)/{2(1+ρ^x,1/2)}1/2. And the p value of GME can be obtained as(6)pGME=2∑x=01‍∫ΘxzH‍∫0∞‍∫L∞‍ϕxzx,z1/2,zHdzxdz1/2dzH+2∫Θ1/2(z)‍Φ(-t∧0)+ρ1/2z1-ρ1/221/2dΦ(z),where L=t{2(1+ρ^x,1/2)}1/2-z1/2 for t=zGME.


### 2.3.p Value of Replication Data Using the Robust Method

In the discovery stage, thep value of robust association tests, including MAX3, MIN2, ZGMS, and ZGME, can be obtained as described in Section 2.2. For the p 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. LetZx(i) for x=0,1/2,1 be the CATT optimal for recessive, additive, and dominant models and let Px(i) be corresponding p value for ith stage (i=1 for discovery and i=2 for replication stages). Also, denote Zx(i)*=(Zx(i)+Z1/2(i))/{2(1+ρ^x,1/2(i))}1/2 for x=0,1/2,1. Then, for CATT with a preselected genetic model, Px(2)=1-Φ(sign(Zx(1)·Zx(2))·|Zx(2)|) using a one-sided p value given the direction of association from the discovery stage, and Px(2)*=1-Φ(sign(Zx(1)*·Zx(2)*)·|Zx(2)*|). Moreover, denote the test statistics and p values using Pearson’s Chi-square test from the ith stage as Tchi2(i) and Pchi2(i). Further, let HWDTT from the ith stage be ZH(i). Then, the second stage p values, using MAX3, MIN2, ZGMS, and ZGME, denoted as PMAX3(2), PMIN2(2), PGMS(2), and PGME(2), can be obtained as follows: (7)PMAX3(2)=Px*(2),wherex*=argmin⁡x∈0,1/2,1⁡Px(1),PMIN2(2)=P1/2(2)·IP1/21≤Pchi22+Pchi22·IP1/21>Pchi22,PGMS(2)=P0(2)I(ZH(1)>c)+P1(2)I(ZH(1)<-c)+P1/22IZH1≤c,ifZ1/21>0;=P02I(ZH1<-c)+P12IZH1>c+P1/22IZH1≤c,ifZ1/2(1)≤0,PGME(2)=P0(2)*I(ZH(1)>c)+P1(2)*I(ZH(1)<-c)+P1/22*IZH1≤c,ifZ1/2(1)>0;=P0(2)*I(ZH(1)<-c)+P1(2)*I(ZH(1)>c)+P1/22*IZH1≤c,ifZ1/2(1)≤0.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 stagep values, the p values from the two stages are independent under the null hypothesis. This is because, under the null hypothesis, the probability of Z1/2 being positive or negative is simply 1/2, and the probability of the selection of a certain genetic model is also a constant (αH for the recessive and dominant models and 1-2αH for the additive model).


### 2.4. Combined Test Usingp Values and Its Statistical Significance

For a given robust test, we can consider the joint analysis by combiningp values from the discovery and replication stages of GWAS. We consider using p 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 p value for each data set might be relatively simple, and the distribution of p 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 p values after -2log⁡ transformation; that is, ZFC=-2w1log⁡(P(1))-2w2log⁡(P(2)), where P(i) is p value from i=0 for discovery and i=1 for replication stages using a given robust test. A specification of w1=w2=1 gives the same weight for discovery and replication stages, and one can consider w1=2πs and w2=2(1-πs) where πs=ND/(ND+NR), and ND and NR are sample sizes of the discovery and replication data sets. A linear combination of two P values after taking the inverse of the standard normal cumulative distribution is given by ZLC={w1Φ-1(1-P(1)/2)+w2Φ-1(1-P(2))}/w12+w22 with a natural choice of w1=πs and w2=1-πs. Let the significance level of the discovery stage be αD, which means that markers with P(1)<αD are selected and replicated in the replication stage. The p value of combined test can then be obtained as pFC=PH0(P(1)<αD,ZFC>zFC) where the observed value of ZFC is zFC. The pFC are calculated as e-zFC/2(1+zFC/2+log⁡αD) for equal weights where zFC>-2log⁡αD and (w1/w1-w2)e-zFC/2w1-(w2/(w1-w2))e-zFC/2w2αd-(w1-w2)/w2 for unequal weights where zFC>-2w1log⁡αD. Detailed derivations are described in the Appendix. Equivalently, for an overall type I error threshold for a single marker of α, one may obtain the threshold CFC of ZFC that satisfies PH0(P(1)<αD,ZFC>CFC)≤α. Similarly, for the ZLC, the p value is calculated as pLC=PH0(P1<αD,ZLC>zLC) = ∫z1-αD/2∞‍ϕ(z)[1-Φ((w12+w22zLC-w1z)/w2)]dz for zLC>z1-αD/2 where the observed value of ZLC=zLC.


## 2.1. Data and Notation

For a marker with two allelesA and B, let the frequencies of B in cases and controls be p=P(B∣case) and q=P(B∣control). Denote three genotypes by G0=AA, G1=AB, and G2=BB. In case-control association studies, r cases and s controls are independently sampled from each population. The observed genotype counts for (G0,G1,G2) are (r0,r1,r2) in the cases and (s0,s1,s2) in the controls. Disease prevalence is denoted by k=P(disease) and penetrance by fi=P(disease∣Gi) for i=0,1,2. Two genotype relative risks (GRRs) are denoted by λ1=f1/f0 and λ2=f2/f0 using f0>0 as baseline penetrance. Under the null hypothesis of no association H0:f0=f1=f2=k or alternatively H0:λ2=λ1=1. Genetic model is recessive (REC), additive (ADD), multiplicative (MUL), and dominant (DOM) when λ1=1, λ1=(1+λ2)/2, λ1=λ21/2, and λ2=λ1, 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 underH0. The test is given by (1)Tchi2=∑j=02‍rj-njr/n2njr/n+∑j=02‍sj-njs/n2njs/n,where ni=ri+si for i=0,1,2 and n=r+s. Under Hardy-Weinberg equilibrium (HWE), an allele-based test (ABT) and CATT with scores (0,x,1) for (G0,G1,G2), where 0≤x≤1, are given by (2)ZABT=n1/22r2s0+s1-2s2r0+r12rs2n0+n1n1+2n21/2,Zx=n1/2∑i=02xisri-rsirsnn∑i=02xi2ni-∑i=02xini21/2,where (x0,x1,x2)=(0,x,1) [9]. The optimal choices of x for the recessive (REC), additive/multiplicative (ADD/MUL), and dominant (DOM) models are x=0,1/2 and 1, respectively [9, 20]. Both Zx and ZABT asymptotically follow a standard normal distribution under H0. Zx can be used even when HWE does not hold. However, without the HWE assumption, ZABT 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 MAX3=max⁡(|Z0|,|Z1/2|,|Z1|). Parametric bootstrap or permutation methods can be used to find the p value of MAX3 [4].Let thep values of Pearson’s Chi-square test and CATT under the additive genetic model Z1/2 be Pchi2 and P1/2, respectively. WTCCC [5] proposed an alternative robust test MIN2=min⁡(Pchi2,P1/2). Joo et al. [15] derived the asymptotic null distribution of MIN2 and using their result the p value of MIN2 can be obtained as (3)PMIN2=12exp⁡-12H1-1(1-MIN2)+12MIN2-12π∫H1-11-MIN2-2log⁡⁡MIN2‍e-v/2arcsin2H1-11-MIN2v-1dv,where H1 and H2 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 (4)ZH=rs/n1/2(Δ^p-Δ^q)1-n2/n-n1/(2n)n2/n+n1/(2n),where Δ^p=p^2-(p^2+p^1/2)2 and Δ^q=q^2-(q^2+q^1/2)2 are the estimates of Δp and Δq, where p^i=ri/r and q^i=si/s. Here, Δ denotes the Hardy-Weinberg disequilibrium (HWD) coefficient defined by Pr(BB)-{Pr(AB)/2+Pr(BB)}2 and Δp and Δq denote the HWD coefficient in cases and controls, respectively.Zheng and Ng [13] used the information contained in the signs of (Δp,Δq) to determine the genetic models in their two-phase method. Their two-phase statistic ZGMS is given by ZGMS=Z0 if ZH>c, Z1 if ZH<-c, and Z1/2 otherwise, where c=Φ-1(1-αH) for αH=0.05. The asymptotic correlations between ZH 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 B is the risk allele, Joo et al. [14] studied the behavior of ZGMS when either one of the alleles can be a risk allele. They chose the risk allele based on the sign of Z1/2; that is, if Z1/2>0, B is the risk allele, and Z0, Z1/2, and Z1 are chosen for REC, ADD, and DOM models, respectively. If Z1/2<0, the respective test statistics are chosen to be -Z1, -Z1/2, and -Z0. They incorporate this property in defining the test statistic for genetic model selection (ZGMS) and calculating the p value. Let Θ0(z)={z:z>c}, Θ1/2(z)={z:|z|<c}, and Θ1(z)={z:z<-c}. Then, the p value of this method can be obtained by(5)pGMS=2∑x=01‍∫Θx(zH)‍∫0∞‍∫t∞‍ϕx(zx,z1/2,zH)dzxdz1/2dzH+2∫Θ1/2(z)‍Φ(-t∧0)+ρ1/2z1-ρ1/221/2dΦ(z),where ρx=Corr(Zx,ZH) in (5) and ρx,1/2=Corr(Zx,Z1/2) (x=0,1) are replaced by their consistent estimates. Here, t=zGMS and (-t∧0)=min⁡(-t,0). Moreover, zGMS and z1/2 are the observed values of ZGMS and Z1/2, respectively.While studying the properties of GMS, Joo et al. [14] noticed that the probability of selecting the true recessive or dominant models using ZH 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 ZGME which is the same as the ZGMS described above except Zx for x=0,1/2,1 is replaced by Zx* where Zx*=(Zx+Z1/2)/{2(1+ρ^x,1/2)}1/2. And the p value of GME can be obtained as(6)pGME=2∑x=01‍∫ΘxzH‍∫0∞‍∫L∞‍ϕxzx,z1/2,zHdzxdz1/2dzH+2∫Θ1/2(z)‍Φ(-t∧0)+ρ1/2z1-ρ1/221/2dΦ(z),where L=t{2(1+ρ^x,1/2)}1/2-z1/2 for t=zGME.


## 2.3.p Value of Replication Data Using the Robust Method

In the discovery stage, thep value of robust association tests, including MAX3, MIN2, ZGMS, and ZGME, can be obtained as described in Section 2.2. For the p 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. LetZx(i) for x=0,1/2,1 be the CATT optimal for recessive, additive, and dominant models and let Px(i) be corresponding p value for ith stage (i=1 for discovery and i=2 for replication stages). Also, denote Zx(i)*=(Zx(i)+Z1/2(i))/{2(1+ρ^x,1/2(i))}1/2 for x=0,1/2,1. Then, for CATT with a preselected genetic model, Px(2)=1-Φ(sign(Zx(1)·Zx(2))·|Zx(2)|) using a one-sided p value given the direction of association from the discovery stage, and Px(2)*=1-Φ(sign(Zx(1)*·Zx(2)*)·|Zx(2)*|). Moreover, denote the test statistics and p values using Pearson’s Chi-square test from the ith stage as Tchi2(i) and Pchi2(i). Further, let HWDTT from the ith stage be ZH(i). Then, the second stage p values, using MAX3, MIN2, ZGMS, and ZGME, denoted as PMAX3(2), PMIN2(2), PGMS(2), and PGME(2), can be obtained as follows: (7)PMAX3(2)=Px*(2),wherex*=argmin⁡x∈0,1/2,1⁡Px(1),PMIN2(2)=P1/2(2)·IP1/21≤Pchi22+Pchi22·IP1/21>Pchi22,PGMS(2)=P0(2)I(ZH(1)>c)+P1(2)I(ZH(1)<-c)+P1/22IZH1≤c,ifZ1/21>0;=P02I(ZH1<-c)+P12IZH1>c+P1/22IZH1≤c,ifZ1/2(1)≤0,PGME(2)=P0(2)*I(ZH(1)>c)+P1(2)*I(ZH(1)<-c)+P1/22*IZH1≤c,ifZ1/2(1)>0;=P0(2)*I(ZH(1)<-c)+P1(2)*I(ZH(1)>c)+P1/22*IZH1≤c,ifZ1/2(1)≤0.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 stagep values, the p values from the two stages are independent under the null hypothesis. This is because, under the null hypothesis, the probability of Z1/2 being positive or negative is simply 1/2, and the probability of the selection of a certain genetic model is also a constant (αH for the recessive and dominant models and 1-2αH for the additive model).


## 2.4. Combined Test Usingp Values and Its Statistical Significance

For a given robust test, we can consider the joint analysis by combiningp values from the discovery and replication stages of GWAS. We consider using p 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 p value for each data set might be relatively simple, and the distribution of p 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 p values after -2log⁡ transformation; that is, ZFC=-2w1log⁡(P(1))-2w2log⁡(P(2)), where P(i) is p value from i=0 for discovery and i=1 for replication stages using a given robust test. A specification of w1=w2=1 gives the same weight for discovery and replication stages, and one can consider w1=2πs and w2=2(1-πs) where πs=ND/(ND+NR), and ND and NR are sample sizes of the discovery and replication data sets. A linear combination of two P values after taking the inverse of the standard normal cumulative distribution is given by ZLC={w1Φ-1(1-P(1)/2)+w2Φ-1(1-P(2))}/w12+w22 with a natural choice of w1=πs and w2=1-πs. Let the significance level of the discovery stage be αD, which means that markers with P(1)<αD are selected and replicated in the replication stage. The p value of combined test can then be obtained as pFC=PH0(P(1)<αD,ZFC>zFC) where the observed value of ZFC is zFC. The pFC are calculated as e-zFC/2(1+zFC/2+log⁡αD) for equal weights where zFC>-2log⁡αD and (w1/w1-w2)e-zFC/2w1-(w2/(w1-w2))e-zFC/2w2αd-(w1-w2)/w2 for unequal weights where zFC>-2w1log⁡αD. Detailed derivations are described in the Appendix. Equivalently, for an overall type I error threshold for a single marker of α, one may obtain the threshold CFC of ZFC that satisfies PH0(P(1)<αD,ZFC>CFC)≤α. Similarly, for the ZLC, the p value is calculated as pLC=PH0(P1<αD,ZLC>zLC) = ∫z1-αD/2∞‍ϕ(z)[1-Φ((w12+w22zLC-w1z)/w2)]dz for zLC>z1-αD/2 where the observed value of ZLC=zLC.


## 3. Simulation Results

### 3.1. Type I Error

Table1 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 (πs) of 0.5 and 0.6 were considered for the minor allele frequency (MAF) of 0.3 and 0.4. We considered M=10 markers to control the genome-wide false positive rate at α=0.05 with the Bonferroni correction. We did not consider a larger M 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 M=10, we performed 20,000 simulations which result in less than 10% of a coefficient of variation for a significance level 0.05/M=0.005 for each marker [24]. The test statistics considered are Z1/2, Pearson’s Chi-square test, MIN2, MAX3, GMS, and GME. For the second stage analysis, we considered a replication-based analysis, ZFC, and ZLC as proposed above. The results are based on the situation under HWE (HWE coefficient F=0). 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 F=0.05 (results not shown).Table 1
Type I error rates of three approaches—replication-based (REP) test, Fisher’s combination (ZFC), and linear combination of test (ZLC)—based on the CATT with an additive model (Z1/2), χ2, MAX3, MIN2, GMS, and GME. The disease prevalence K=0.1, M=10 markers, r=1,500 cases, and s=1,500 controls are considered based on 20,000 simulations.
MAF
π
s
α
D
F
=
0
Z
1
/
2
χ
2
MAX3
MIN2
GMS
GME
0.3
0.5
0.05
REP
0.00530
0.00455
0.00505
0.00485
0.0050
0.00490
Z
FC
0.00500
0.00535
0.00495
0.00510
0.00515
0.00460
Z
LC
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
Z
FC
0.00565
0.00535
0.00565
0.00545
0.00565
0.00540
Z
LC
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
Z
FC
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
Z
FC
0.00465
0.00490
0.00485
0.00500
0.00455
0.00460
Z
LC
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
Z
FC
0.00575
0.00430
0.00460
0.00535
0.00460
0.00500
Z
LC
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
Z
FC
0.00515
0.00510
0.00495
0.00475
0.00540
0.00500
Z
LC
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
Z
FC
0.00460
0.00470
0.00420
0.00490
0.00455
0.00440
Z
LC
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
Z
FC
0.00520
0.00370
0.00495
0.00450
0.00515
0.00510
Z
LC
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 Figure1, we considered M=10 markers, a disease prevalence of 10%, the same genotype relative risk for two stages (r1=1.4 and r2=1.4), and 1,000 cases and 1,000 controls. 2,000 simulations were performed under HWE (F=0) to control the genome-wide false positive rate at α=0.05. 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 r1=1.6 and r2=1.4, 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 Z1/2, χ2, and MIN2 are less powerful. Under the additive model, Z1/2 is most powerful, as expected, and χ2 is least powerful. Other robust methods perform well with a slight decrease in power compared to Z1/2. 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 forM=10 markers, genotype relative risks of both stages = 1.4, and disease prevalence K=0.1 under the recessive, additive, and dominant models. 1,000 cases and 1,000 controls are considered to control α=0.05. The first stage type I error rate for discovery is αD=0.05. Six test statistics, Z1/2, χ2, MAX3, MIN2, GMS, and GME, are considered. The first, second, and third columns depict powers using the replication-based test, ZFC, and ZLC, respectively.Figure 2
Empirical powers based on 2,000 simulations forM=10 markers; genotype relative risks of two stages are different (r1=1.6, r2=1.4); disease prevalence K=0.1 under the recessive, additive, and dominant models. 1,000 cases and 1,000 controls are considered to control α=0.05. The first stage type I error rate for discovery is αD=0.05. Six test statistics, Z1/2, χ2, MAX3, MIN2, GMS, and GME, are considered. The first, second, and third columns depict powers using the replication-based test, ZFC, and ZLC, respectively.


## 3.1. Type I Error

Table1 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 (πs) of 0.5 and 0.6 were considered for the minor allele frequency (MAF) of 0.3 and 0.4. We considered M=10 markers to control the genome-wide false positive rate at α=0.05 with the Bonferroni correction. We did not consider a larger M 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 M=10, we performed 20,000 simulations which result in less than 10% of a coefficient of variation for a significance level 0.05/M=0.005 for each marker [24]. The test statistics considered are Z1/2, Pearson’s Chi-square test, MIN2, MAX3, GMS, and GME. For the second stage analysis, we considered a replication-based analysis, ZFC, and ZLC as proposed above. The results are based on the situation under HWE (HWE coefficient F=0). 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 F=0.05 (results not shown).Table 1
Type I error rates of three approaches—replication-based (REP) test, Fisher’s combination (ZFC), and linear combination of test (ZLC)—based on the CATT with an additive model (Z1/2), χ2, MAX3, MIN2, GMS, and GME. The disease prevalence K=0.1, M=10 markers, r=1,500 cases, and s=1,500 controls are considered based on 20,000 simulations.
MAF
π
s
α
D
F
=
0
Z
1
/
2
χ
2
MAX3
MIN2
GMS
GME
0.3
0.5
0.05
REP
0.00530
0.00455
0.00505
0.00485
0.0050
0.00490
Z
FC
0.00500
0.00535
0.00495
0.00510
0.00515
0.00460
Z
LC
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
Z
FC
0.00565
0.00535
0.00565
0.00545
0.00565
0.00540
Z
LC
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
Z
FC
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
Z
FC
0.00465
0.00490
0.00485
0.00500
0.00455
0.00460
Z
LC
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
Z
FC
0.00575
0.00430
0.00460
0.00535
0.00460
0.00500
Z
LC
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
Z
FC
0.00515
0.00510
0.00495
0.00475
0.00540
0.00500
Z
LC
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
Z
FC
0.00460
0.00470
0.00420
0.00490
0.00455
0.00440
Z
LC
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
Z
FC
0.00520
0.00370
0.00495
0.00450
0.00515
0.00510
Z
LC
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 Figure1, we considered M=10 markers, a disease prevalence of 10%, the same genotype relative risk for two stages (r1=1.4 and r2=1.4), and 1,000 cases and 1,000 controls. 2,000 simulations were performed under HWE (F=0) to control the genome-wide false positive rate at α=0.05. 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 r1=1.6 and r2=1.4, 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 Z1/2, χ2, and MIN2 are less powerful. Under the additive model, Z1/2 is most powerful, as expected, and χ2 is least powerful. Other robust methods perform well with a slight decrease in power compared to Z1/2. 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 forM=10 markers, genotype relative risks of both stages = 1.4, and disease prevalence K=0.1 under the recessive, additive, and dominant models. 1,000 cases and 1,000 controls are considered to control α=0.05. The first stage type I error rate for discovery is αD=0.05. Six test statistics, Z1/2, χ2, MAX3, MIN2, GMS, and GME, are considered. The first, second, and third columns depict powers using the replication-based test, ZFC, and ZLC, respectively.Figure 2
Empirical powers based on 2,000 simulations forM=10 markers; genotype relative risks of two stages are different (r1=1.6, r2=1.4); disease prevalence K=0.1 under the recessive, additive, and dominant models. 1,000 cases and 1,000 controls are considered to control α=0.05. The first stage type I error rate for discovery is αD=0.05. Six test statistics, Z1/2, χ2, MAX3, MIN2, GMS, and GME, are considered. The first, second, and third columns depict powers using the replication-based test, ZFC, and ZLC, respectively.


## 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 Z1/2. In the replication stage, 168 SNPs with p value less than 1×10-4 in the first stage based on Z1/2 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 Z1/2 since MIN2 produces stronger evidence for the additional SNPs than Z1/2 does. The Manhattan plots of using MIN2 and Z1/2 are presented in Figure 3. One example is rs385272 located in chromosome 2, which had a p value of 1.37×10-7 which reached significance level at Bonferroni correction in discovery samples alone, whereas Z1/2 yielded a p value greater than 1×10-4. 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 Z1/2 (b). The x axis is chromosomal location and the y axis is the significance (-log⁡10P) of association. The horizontal line corresponds to the significance level 10-4.
(a)
(b)Since we do not have replication data for these additional SNPs selected using MIN2 because the first stage selection was based onZ1/2 in Yoon et al. [25], just for illustration purpose of the proposed methods, we present the results of three SNPs including rs2131877 that was reported by Yoon et al. [25]. When the significance level in the discovery stage is set at αD=5×10-5 so that all these exemplary SNPs can be selected in the discovery stage; the p value of combined test based on four robust methods (MAX3, MIN2, GMS, and GME) as well as Z1/2 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 rs2131877 was found to be significant with Bonferroni correction (p value <2.03×10-7) by all except MAX3 method.Table 2
For selected exemplary three SNPs for testing association with NSCLC,p value of combined test using additive CATT (Z1/2), Pearson’s Chi-square test (Tchi2), MAX3, MIN2, ZGMS, and ZGME.
SNP
p value of Z1/2
p value of Tchi2
Discovery
Replication
Combined test
Discovery
Replication
Combined test
rs2131877
7.88
×
1
0
-
5
1.04
×
1
0
-
4
7.97
×
1
0
-
8
1.40
×
1
0
-
4
1.49
×
1
0
-
4
1.84
×
1
0
-
7
rs905551
1.83
×
1
0
-
5
7.02
×
1
0
-
3
7.70
×
1
0
-
6
8.06
×
1
0
-
5
4.89
×
1
0
-
2
1.40
×
1
0
-
5
rs1695109
2.48
×
1
0
-
4
3.46
×
1
0
-
2
2.17
×
1
0
-
6
4.56
×
1
0
-
5
1.53
×
1
0
-
1
2.07
×
1
0
-
5
SNP
p value of MAX3
p value of MIN2
Discovery
Replication
Combined test
Discovery
Replication
Combined test
rs2131877
1.53
×
1
0
-
4
4.05
×
1
0
-
2
1.92
×
1
0
-
5
1.32
×
1
0
-
4
1.04
×
1
0
-
4
1.26
×
1
0
-
7
rs905551
4.50
×
1
0
-
5
7.02
×
1
0
-
3
1.92
×
1
0
-
6
1.34
×
1
0
-
4
4.89
×
1
0
-
2
1.99
×
1
0
-
5
rs1695109
3.54
×
1
0
-
5
2.63
×
1
0
-
2
4.64
×
1
0
-
6
2.36
×
1
0
-
5
2.63
×
1
0
-
2
3.35
×
1
0
-
6
SNP
p value of ZGMS
p value of ZGME
Discovery
Replication
Combined test
Discovery
Replication
Combined test
rs2131877
1.86
×
1
0
-
4
1.04
×
1
0
-
4
1.71
×
1
0
-
7
1.03
×
1
0
-
4
1.04
×
1
0
-
4
1.02
×
1
0
-
7
rs905551
5.19
×
1
0
-
5
7.02
×
1
0
-
3
2.16
×
1
0
-
6
7.35
×
1
0
-
5
8.01
×
1
0
-
3
3.20
×
1
0
-
6
rs1695109
6.89
×
1
0
-
4
1.27
×
1
0
-
1
3.85
×
1
0
-
5
2.69
×
1
0
-
5
4.19
×
1
0
-
2
5.40
×
1
0
-
6


## 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 byp 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 thep value of each stage might be relatively simple. Because the asymptotic distribution of the p value under the null hypothesis of no association is easy to handle, the combined test using p 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 p 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 p 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 p 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 p 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.



---

*Source: 461593-2015-08-04.xml*