# Robust Joint Analysis with Data Fusion in Two-Stage Quantitative Trait Genome-Wide Association Studies

Authors: Dong-Dong Pan; Wen-Jun Xiong; Ji-Yuan Zhou; Ying Pan; Guo-Li Zhou; Wing-Kam Fung
Journal: Computational and Mathematical Methods in Medicine
 (2013)

---

## Abstract

Genome-wide association studies (GWASs) in identifying the disease-associated genetic variants
have been proved to be a great pioneering work. Two-stage design and analysis are often adopted in
GWASs. Considering the genetic model uncertainty, many robust procedures have been proposed and
applied in GWASs. However, the existing approaches mostly focused on binary traits, and few work
has been done on continuous (quantitative) traits, since the statistical significance of these robust tests
is difficult to calculate. In this paper, we develop a powerfulF-statistic-based robust joint analysis
method for quantitative traits using the combined raw data from both stages in the framework of
two-staged GWASs. Explicit expressions are obtained to calculate the statistical significance and
power. We show using simulations that the proposed method is substantially more robust than the
F-test based on the additive model when the underlying genetic model is unknown. An example for
rheumatic arthritis (RA) is used for illustration.


---

## Body

## 1. Introduction

Genome-wide association studies (GWASs) have identified a large number of genomic regions (especially single-nucleotide polymorphisms (SNPs)) with a wide variety of complex traits/diseases. In a GWAS, two most common types of data, qualitative (or binary) and quantitative (or continuous) traits, are analyzed and two contentious points are often faced; one is how to construct the test statistic considering the genetic model uncertainty and the other is how to evaluate the statistical significance for controlling the false positive rates efficiently (e.g., [1, 2]). Considering these issues, a lot of work has been done on the binary trait in the past 10 years (e.g., [3–7]). Computer algorithms have also been developed to calculated the significance level of robust tests in GWASs, taking into account the genetic model uncertainty [8]. However, few work has been done on continuous traits, only recently So and Sham [9] proposed a MAX3 based on score test statistics, and Li et al. [10] gave a MAX3 based on F-test statistics. Note that these tests just focus on single-marker analysis in one-stage analysis.Although the costs of whole-genome genotyping are decreasing with the high-throughput biological technology, the total costs for a GWAS are still very expensive due to the thousands of sampling units and huge amounts of single-nucleotide polymorphisms. In order to save the costs, the two-stage design and the corresponding statistical analysis where all the SNPs are genotyped in Stage 1 on a portion of the samples and the promising SNPs with smallP-values (e.g., <0.001) based on some efficient tests are further screened on the remaining subjects, are often adopted in practice (e.g., [11–15]).In genetic association studies, especially GWASs, genetic markers are routinely tested under the assumption of additive effects. Although convenient to use, those tests are optimal only when the true underlying genetic model is additive so that they are not robust against the genetic model misspecification. To our best knowledge, few work has been done on the two-stage joint analysis for quantitative trait GWASs allowing for genetic model uncertainty. Here, we attempt to develop a joint analysis method with data fusion in the two-stage design usingF-statistic, since F-test is commonly employed from the linear regression model for quantitative trait, and Li et al. [10] show that MAX3 based on F-statistics is more powerful than So and Sham’s method by extensively numerical simulation.The content of this paper is organized as follows. In Section2, we give some notations and the proposed robust joint test statistics. Further, we derive the asymptotic distribution of the test statistics under the null and the alternative hypotheses. In Section 3, we show that the proposed joint analysis method is substantially more robust than the additive-model-based F-test from the numerical results of power comparison when the real genetic model is unknown. After that, an illustrative example for rheumatic arthritis (RA) is presented. Finally, we give some discussion of this paper in Section 4.


## 2. Methods

### 2.1. Notations

Assume thatn individuals are randomly selected to be genotyped in a two-staged GWAS for a certain quantitative trait and that π is the sampling proportion in Stage 1. Let n1=nπ and n2=n(1-π) be the sample sizes for Stages 1 and 2, respectively. Consider a biallelic marker with two alleles G and g. Without loss of generality, we assume that G is the minor or high-risk allele. We suppose that the total m SNPs are genotyped on the samples of Stage 1, and SNPs with P-values less than γ in Stage 1 will be further genotyped and tested in Stage 2. Let the significance level be α, and then the genome-wide significance level per SNP is α/m with the Bonferroni adjustments. Let Y1=(y1,y2,…,yn1)′ and Y2=(yn1+1,yn1+2,…,yn)′ be the observed quantitative outcome vectors for Stage 1 and Stage 2, respectively. Without loss of generality, we assume that the first n10 individuals in Stage 1 have the genotype gg, the second n11 individuals in Stage 1 have the genotype Gg, and the last n12 subjects in Stage 1 possess the genotype GG. Similarly, the first n20 subjects in Stage 2 have the genotype gg, the second n21 individuals in Stage 2 have the genotype Gg, and the last n22 subjects in Stage 2 possess the genotype GG. Let 0k=(0,0,…,0)k×1′ and 1k=(1,1,…,1)k×1′, and let Ok×j be the k×j matrix with all its entries being zero and In be the n×n identity matrix.


### 2.2.F-Statistic-Based Robust Joint Analysis

We firstly briefly introduceF-statistic-based MAX3 by Li et al. [10] just using the data from Stage 1. Consider the following linear regression model:
(1)yi=β0+giβ1+εi,εi~N(0,σ2),i=1,2,…,n1,
where β0 is the nuisance parameter for the intercept, β1 is the parameter of interest for genetic effect, and gi is the genotype value, which takes 0, 1, or 2 corresponding to the count of G at a marker locus for the ith subject, i=1,2,…,n1. The hypotheses of interest are
(2)H0:β1=0⟷H1:β1≠0.
The variable gi in the previously stated equation is coded differently for the three common genetic models. Let X1R=(1n1,G1R), X1A=(1n1,G1A), and X1D=(1n1,G1D) be the design matrices under three commonly used genetic models, where G1R=(0n10+n11′,1n12′)′ corresponds to the recessive model, G1A=(g1,g2,…,gn1)′ corresponds to the additive model, and G1D=(0n10′,1n11+n12′)′ is for the dominant model. Denote X1=(1n1,x11,x12), where x11=(0n10′,1n11′,0n12′)′ and x12=(0n10′,0n11′,1n12′)′. The modified F-test statistics under the recessive, additive, and dominant models for Stage 1 are given by
(3)F1R=Y1′[X1R(X1R′X1R)-1X1R′-1n1(1n1′1n1)-11n1′]Y1Y1′[In1-X1(X1′X1)-1X1′]Y1/(n1-3)=(Z1R)2RSS1/(n1-3),F1A=Y1′[X1A(X1A′X1A)-1X1A′-1n1(1n1′1n1)-11n1′]Y1Y1′[In1-X1(X1′X1)-1X1′]Y1/(n1-3)=(Z1A)2RSS1/(n1-3),F1D=Y1′[X1D(X1D′X1D)-1X1D′-1n1(1n1′1n1)-11n1′]Y1Y1′[In1-X1(X1′X1)-1X1′]Y1/(n1-3)=(Z1D)2RSS1/(n1-3),
where
(4)Z1R=(n10+n11)n12n1(y-n10+n11-y-n12),Z1D=n10(n11+n12)n1(y-n10-y-n11+n12),Z1A=(n10(n11+2n12)y-n10-n11(n10-n12)y-n11-n12(2n10+n11)y-n12)×(n1[n10(n11+4n12)+n11n12])-1,y-n10=1n10∑j=1n10‍yj,y-n11=1n11∑j=n10+1n10+n11‍yj,y-n12=1n12∑j=n10+n11+1n1‍yj,y-n10+n11=1n10+n11∑j=1n10+n11‍yj,y-n11+n12=1n11+n12∑j=n10+1n1‍yj.
The robust test statistic in Stage 1 is
(5)F1MAX=max{F1R,F1A,F1D}.We now give the proposed robust joint analysis. In the framework of two-stage design GWAS of quantitative traits, the SNPs withP-values less than γ will be genotyped on the remaining n2 subjects in Stage 2. Following the previous notation for Stage 1, corresponding to the recessive, additive, and dominant models, the genotype data in Stage 2 are denoted by G2R=(0n20+n21′,1n22′)′, G2A=(gn1+1,gn1+2,…,gn)′, and G2D=(0n20′,1n21+n22′)′, respectively, and the design matrices are X2R=(1n2,G2R), X2A=(1n2,G2A), and X2D=(1n2,G2D), respectively. Denote X2=(1n2,x21,x22), where x21=(0n20′,1n21′,0n22′)′ and x22=(0n20′,0n21′,1n22′)′. Then, we can obtain three modified F-test statistics under the recessive, additive, and dominant models for Stage 2 similarly, and denote them by F2R, F2A, and F2D. Let Y=(Y1′,Y2′)′, GR=(G1R′,G2R′)′, GA=(G1A′,G2A′)′, and GD=(G1D′,G2D′)′. Denote N0=n10+n20, N1=n11+n21, and N2=n12+n22 for the combined sample sizes from two stages, corresponding to three genotypes. Then the proposed F-test statistics under three genetic models on the basis of the combined data are as follows:
(6)FJR=Y′[XR(XR′XR)-1XR′-1n(1n′1n)-11n′]YY′[In-W(W′W)-1W′]Y/(n-6)=(ZJR)2RSSJ/(n-6),FJA=Y′[XA(XA′XA)-1XA′-1n(1n′1n)-11n′]YY′[In-W(W′W)-1W′]Y/(n-6)=(ZJA)2RSSJ/(n-6),FJD=Y′[XD(XD′XD)-1XD′-1n(1n′1n)-11n′]YY′[In-W(W′W)-1W′]Y/(n-6)=(ZJD)2RSSJ/(n-6),
where XR=(1n,GR), XA=(1n,GA), XD=(1n,GD), and   W=(X1On1×3On2×3X2),
(7)ZJR=(N0+N1)N2n(y-01-y-2),ZJD=N0(N1+N2)n(y-0-y-12),ZJA=(N0(N1+2N2)y-0-N1(N0-N2)y-1-N2(2N0+N1)y-2)×(n[N0(N1+4N2)+N1N2])-1,y-0=1N0(∑j=1n10‍yj+∑j=n1+1n1+n20‍yj),y-1=1N1(∑j=n10+1n10+n11‍yj+∑j=n1+n20+1n1+n20+n21‍yj),y-2=1N2(∑j=n10+n11+1n1‍yj+∑j=n1+n20+n21+1n‍yj),y-01=N0y-0+N1y-1N0+N1,y-12=N1y-1+N2y-2N1+N2.
Furthermore, we propose the joint testing statistic as
(8)FJMAX=max{FJR,FJA,FJD}.In order to calculate the power of the proposed joint analysis, we have to get the thresholds, which is determined by the significance level. Denote the threshold for choosing the promising SNPs in Stage 1 byu1, which is the solution of
(9)PrH0(F1MAX>u1)=γ.
Since the genome-wide significance level is α/m, in order to control the false positive rate, we have
(10)PrH0(F1MAX>u1,FJMAX>uJ)=α/m,
where uJ is the cut-off point for the joint statistic. Once we have u1 and uJ, the power is calculated by
(11)PrH1(F1MAX>u1,FJMAX>uJ).We now give the detail to calculate the cut-off point and power above. The left side of (10) can be further expressed as
(12)PrH0(F1MAX>u1,FJMAX>uJ)=1-PrH0(F1MAX≤u1)-PrH0(FJMAX≤uJ)+PrH0(F1MAX≤u1,FJMAX≤uJ).
For controlling the type I error rate and calculating the power, we need to know the distribution or the asymptotic distribution of (Z1R,Z1A,Z1D,ZJR,ZJA,ZJD,RSS1,RRSJ)′ under both H0 and H1.Note that whetherH0 or H1 holds, RSS1 and RSSJ and (Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′ are mutually independent (the proof is given in Appendix A). Denote the correlation matrix of (Z1R,Z1A,Z1D)′ by V1=(vkl)3×3, whose entries are v11=v22=v33=1, v12=v21=n12(2n10+n11)/(n10+n11)[n10(n11+4n12)+n11n12], v13=v31=n10n12/(n10+n11)(n11+n12), and v23=v32=n10(2n12+n11)/(n11+n12)[n10(n11+4n12)+n11n12], respectively. Similarly, let VJ=(vkl*)3×3 be the correlation matrix of (ZJR,ZJA,ZJD)′ with v11*=v22*=v33*=1, v12*=v21*=N2(2N0+N1)/(N0+N1)[N0(N1+4N2)+N1N2], v13*=v31*=N0N2/(N0+N1)(N1+N2), and v23*=v32*=N0(2N2+N1)/(N1+N2)[N0(N1+4N2)+N1N2]. Then, we can derive that RRS1/σ2~χn1-32, RSSJ/σ2~χn-62, and
(13)(Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′|H0~N6(06,σ2(V1ρρ′VJ)),
where ρ=(ρkl)3×3 is the correlation matrix between (Z1R,Z1A,Z1D)′ and (ZJR,ZJA,ZJD)′, with
(14)ρ11=Corr(Z1R,ZJR)=n(n10+n11)n12n1(N0+N1)N2,ρ12=Corr(Z1R,ZJA)=nn12(2n10+n11)n1(n10+n11)[N0(N1+4N2)+N1N2],ρ13=Corr(Z1R,ZJD)=nn12n10n1(n10+n11)N0(N1+N2),ρ21=Corr(Z1A,ZJR)=nn12(2n10+n11)n1[n10(n11+4n12)+n11n12](N0+N1)N2,ρ22=Corr(Z1A,ZJA)=n[n10(n11+4n12)+n11n12]n1[N0(N1+4N2)+N1N2],ρ23=Corr(Z1A,ZJD)=nn10(n11+2n12)n1[n10(n11+4n12)+n11n12]N0(N1+N2),ρ31=Corr(Z1D,ZJR)=nn10n12n1(n11+n12)(N0+N1)N2,ρ32=Corr(Z1D,ZJA)=nn10(n11+2n12)n1(n11+n12)[N0(N1+4N2)+N1N2],ρ33=Corr(Z1D,ZJD)=nn10(n11+n12)n1N0(N1+N2).UnderH1, for a given odds ratio OR=exp(β1) for subjects with two copies of risk allele corresponding to recessive model or one copy of risk allele corresponding to additive or dominant models, we have the following:(i)
when the true genetic model is recessive,(15)(Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′|H1~N6(μR,σ2(V1ρρ′VJ)),whereμR=(μ1RR,μ1RA,μ1RD,μJRR,μJRA,μJRD)′ with
(16)μ1RR=-(n10+n11)n12n1β1,μ1RA=-n12(2n10+n11)β1n1[n10(n11+4n12)+n11n12],μ1RD=-n10n12β1n1(n11+n12),μJRR=-(N0+N1)N2nβ1,μJRA=-N2(2N0+N1)β1n[N0(N1+4N2)+N1N2],μJRD=-N0N2β1n(N1+N2),(ii)
when the true genetic model is additive,(17)(Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′|H1~N6(μA,σ2(V1ρρ′VJ)),whereμA=(μ1AR,μ1AA,μ1AD,μJAR,μJAA,μJAD)′ with
(18)μ1AR=-n12(n11+2n10)β1n1(n10+n11),μ1AA=-n10(n11+4n12)+n11n12n1β1,μ1AD=-n10(n11+2n12)β1n1(n11+n12),μJAR=-N2(N1+2N0)β1n(N0+N1),μJAA=-N0(N1+4N2)+N1N2nβ1,μJAD=-N0(N1+2N2)β1n(N1+N2),(iii)
when the true genetic model is dominant,(19)(Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′|H1~N6(μD,σ2(V1ρρ′VJ)),whereμD=(μ1DR,μ1DA,μ1DD,μJDR,μJDA,μJDD)′ with
(20)μ1DR=-n12n10β1n1(n10+n11),μ1DA=-n10(n11+2n12)β1n1[n10(n11+4n12)+n11n12],μ1DD=-n10(n11+n12)n1β1,μJDR=-N2N0β1n(N0+N1),μJDA=-N0(N1+2N2)β1n[N0(N1+4N2)+N1N2],μJDD=-N0(N1+N2)nβ1.We develop a method for simplifying the calculations ofPrH0(F1MAX≤u1) and PrH0(FJMAX≤uJ) and PrH0(F1MAX≤u1,FJMAX≤uJ). The details are included in Appendix B, and the calculations of PrH1(F1MAX≤u1) and PrH1(FJMAX≤uJ) and PrH1(F1MAX≤u1,FJMAX≤uJ) are essentially similar.


## 2.1. Notations

Assume thatn individuals are randomly selected to be genotyped in a two-staged GWAS for a certain quantitative trait and that π is the sampling proportion in Stage 1. Let n1=nπ and n2=n(1-π) be the sample sizes for Stages 1 and 2, respectively. Consider a biallelic marker with two alleles G and g. Without loss of generality, we assume that G is the minor or high-risk allele. We suppose that the total m SNPs are genotyped on the samples of Stage 1, and SNPs with P-values less than γ in Stage 1 will be further genotyped and tested in Stage 2. Let the significance level be α, and then the genome-wide significance level per SNP is α/m with the Bonferroni adjustments. Let Y1=(y1,y2,…,yn1)′ and Y2=(yn1+1,yn1+2,…,yn)′ be the observed quantitative outcome vectors for Stage 1 and Stage 2, respectively. Without loss of generality, we assume that the first n10 individuals in Stage 1 have the genotype gg, the second n11 individuals in Stage 1 have the genotype Gg, and the last n12 subjects in Stage 1 possess the genotype GG. Similarly, the first n20 subjects in Stage 2 have the genotype gg, the second n21 individuals in Stage 2 have the genotype Gg, and the last n22 subjects in Stage 2 possess the genotype GG. Let 0k=(0,0,…,0)k×1′ and 1k=(1,1,…,1)k×1′, and let Ok×j be the k×j matrix with all its entries being zero and In be the n×n identity matrix.


## 2.2.F-Statistic-Based Robust Joint Analysis

We firstly briefly introduceF-statistic-based MAX3 by Li et al. [10] just using the data from Stage 1. Consider the following linear regression model:
(1)yi=β0+giβ1+εi,εi~N(0,σ2),i=1,2,…,n1,
where β0 is the nuisance parameter for the intercept, β1 is the parameter of interest for genetic effect, and gi is the genotype value, which takes 0, 1, or 2 corresponding to the count of G at a marker locus for the ith subject, i=1,2,…,n1. The hypotheses of interest are
(2)H0:β1=0⟷H1:β1≠0.
The variable gi in the previously stated equation is coded differently for the three common genetic models. Let X1R=(1n1,G1R), X1A=(1n1,G1A), and X1D=(1n1,G1D) be the design matrices under three commonly used genetic models, where G1R=(0n10+n11′,1n12′)′ corresponds to the recessive model, G1A=(g1,g2,…,gn1)′ corresponds to the additive model, and G1D=(0n10′,1n11+n12′)′ is for the dominant model. Denote X1=(1n1,x11,x12), where x11=(0n10′,1n11′,0n12′)′ and x12=(0n10′,0n11′,1n12′)′. The modified F-test statistics under the recessive, additive, and dominant models for Stage 1 are given by
(3)F1R=Y1′[X1R(X1R′X1R)-1X1R′-1n1(1n1′1n1)-11n1′]Y1Y1′[In1-X1(X1′X1)-1X1′]Y1/(n1-3)=(Z1R)2RSS1/(n1-3),F1A=Y1′[X1A(X1A′X1A)-1X1A′-1n1(1n1′1n1)-11n1′]Y1Y1′[In1-X1(X1′X1)-1X1′]Y1/(n1-3)=(Z1A)2RSS1/(n1-3),F1D=Y1′[X1D(X1D′X1D)-1X1D′-1n1(1n1′1n1)-11n1′]Y1Y1′[In1-X1(X1′X1)-1X1′]Y1/(n1-3)=(Z1D)2RSS1/(n1-3),
where
(4)Z1R=(n10+n11)n12n1(y-n10+n11-y-n12),Z1D=n10(n11+n12)n1(y-n10-y-n11+n12),Z1A=(n10(n11+2n12)y-n10-n11(n10-n12)y-n11-n12(2n10+n11)y-n12)×(n1[n10(n11+4n12)+n11n12])-1,y-n10=1n10∑j=1n10‍yj,y-n11=1n11∑j=n10+1n10+n11‍yj,y-n12=1n12∑j=n10+n11+1n1‍yj,y-n10+n11=1n10+n11∑j=1n10+n11‍yj,y-n11+n12=1n11+n12∑j=n10+1n1‍yj.
The robust test statistic in Stage 1 is
(5)F1MAX=max{F1R,F1A,F1D}.We now give the proposed robust joint analysis. In the framework of two-stage design GWAS of quantitative traits, the SNPs withP-values less than γ will be genotyped on the remaining n2 subjects in Stage 2. Following the previous notation for Stage 1, corresponding to the recessive, additive, and dominant models, the genotype data in Stage 2 are denoted by G2R=(0n20+n21′,1n22′)′, G2A=(gn1+1,gn1+2,…,gn)′, and G2D=(0n20′,1n21+n22′)′, respectively, and the design matrices are X2R=(1n2,G2R), X2A=(1n2,G2A), and X2D=(1n2,G2D), respectively. Denote X2=(1n2,x21,x22), where x21=(0n20′,1n21′,0n22′)′ and x22=(0n20′,0n21′,1n22′)′. Then, we can obtain three modified F-test statistics under the recessive, additive, and dominant models for Stage 2 similarly, and denote them by F2R, F2A, and F2D. Let Y=(Y1′,Y2′)′, GR=(G1R′,G2R′)′, GA=(G1A′,G2A′)′, and GD=(G1D′,G2D′)′. Denote N0=n10+n20, N1=n11+n21, and N2=n12+n22 for the combined sample sizes from two stages, corresponding to three genotypes. Then the proposed F-test statistics under three genetic models on the basis of the combined data are as follows:
(6)FJR=Y′[XR(XR′XR)-1XR′-1n(1n′1n)-11n′]YY′[In-W(W′W)-1W′]Y/(n-6)=(ZJR)2RSSJ/(n-6),FJA=Y′[XA(XA′XA)-1XA′-1n(1n′1n)-11n′]YY′[In-W(W′W)-1W′]Y/(n-6)=(ZJA)2RSSJ/(n-6),FJD=Y′[XD(XD′XD)-1XD′-1n(1n′1n)-11n′]YY′[In-W(W′W)-1W′]Y/(n-6)=(ZJD)2RSSJ/(n-6),
where XR=(1n,GR), XA=(1n,GA), XD=(1n,GD), and   W=(X1On1×3On2×3X2),
(7)ZJR=(N0+N1)N2n(y-01-y-2),ZJD=N0(N1+N2)n(y-0-y-12),ZJA=(N0(N1+2N2)y-0-N1(N0-N2)y-1-N2(2N0+N1)y-2)×(n[N0(N1+4N2)+N1N2])-1,y-0=1N0(∑j=1n10‍yj+∑j=n1+1n1+n20‍yj),y-1=1N1(∑j=n10+1n10+n11‍yj+∑j=n1+n20+1n1+n20+n21‍yj),y-2=1N2(∑j=n10+n11+1n1‍yj+∑j=n1+n20+n21+1n‍yj),y-01=N0y-0+N1y-1N0+N1,y-12=N1y-1+N2y-2N1+N2.
Furthermore, we propose the joint testing statistic as
(8)FJMAX=max{FJR,FJA,FJD}.In order to calculate the power of the proposed joint analysis, we have to get the thresholds, which is determined by the significance level. Denote the threshold for choosing the promising SNPs in Stage 1 byu1, which is the solution of
(9)PrH0(F1MAX>u1)=γ.
Since the genome-wide significance level is α/m, in order to control the false positive rate, we have
(10)PrH0(F1MAX>u1,FJMAX>uJ)=α/m,
where uJ is the cut-off point for the joint statistic. Once we have u1 and uJ, the power is calculated by
(11)PrH1(F1MAX>u1,FJMAX>uJ).We now give the detail to calculate the cut-off point and power above. The left side of (10) can be further expressed as
(12)PrH0(F1MAX>u1,FJMAX>uJ)=1-PrH0(F1MAX≤u1)-PrH0(FJMAX≤uJ)+PrH0(F1MAX≤u1,FJMAX≤uJ).
For controlling the type I error rate and calculating the power, we need to know the distribution or the asymptotic distribution of (Z1R,Z1A,Z1D,ZJR,ZJA,ZJD,RSS1,RRSJ)′ under both H0 and H1.Note that whetherH0 or H1 holds, RSS1 and RSSJ and (Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′ are mutually independent (the proof is given in Appendix A). Denote the correlation matrix of (Z1R,Z1A,Z1D)′ by V1=(vkl)3×3, whose entries are v11=v22=v33=1, v12=v21=n12(2n10+n11)/(n10+n11)[n10(n11+4n12)+n11n12], v13=v31=n10n12/(n10+n11)(n11+n12), and v23=v32=n10(2n12+n11)/(n11+n12)[n10(n11+4n12)+n11n12], respectively. Similarly, let VJ=(vkl*)3×3 be the correlation matrix of (ZJR,ZJA,ZJD)′ with v11*=v22*=v33*=1, v12*=v21*=N2(2N0+N1)/(N0+N1)[N0(N1+4N2)+N1N2], v13*=v31*=N0N2/(N0+N1)(N1+N2), and v23*=v32*=N0(2N2+N1)/(N1+N2)[N0(N1+4N2)+N1N2]. Then, we can derive that RRS1/σ2~χn1-32, RSSJ/σ2~χn-62, and
(13)(Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′|H0~N6(06,σ2(V1ρρ′VJ)),
where ρ=(ρkl)3×3 is the correlation matrix between (Z1R,Z1A,Z1D)′ and (ZJR,ZJA,ZJD)′, with
(14)ρ11=Corr(Z1R,ZJR)=n(n10+n11)n12n1(N0+N1)N2,ρ12=Corr(Z1R,ZJA)=nn12(2n10+n11)n1(n10+n11)[N0(N1+4N2)+N1N2],ρ13=Corr(Z1R,ZJD)=nn12n10n1(n10+n11)N0(N1+N2),ρ21=Corr(Z1A,ZJR)=nn12(2n10+n11)n1[n10(n11+4n12)+n11n12](N0+N1)N2,ρ22=Corr(Z1A,ZJA)=n[n10(n11+4n12)+n11n12]n1[N0(N1+4N2)+N1N2],ρ23=Corr(Z1A,ZJD)=nn10(n11+2n12)n1[n10(n11+4n12)+n11n12]N0(N1+N2),ρ31=Corr(Z1D,ZJR)=nn10n12n1(n11+n12)(N0+N1)N2,ρ32=Corr(Z1D,ZJA)=nn10(n11+2n12)n1(n11+n12)[N0(N1+4N2)+N1N2],ρ33=Corr(Z1D,ZJD)=nn10(n11+n12)n1N0(N1+N2).UnderH1, for a given odds ratio OR=exp(β1) for subjects with two copies of risk allele corresponding to recessive model or one copy of risk allele corresponding to additive or dominant models, we have the following:(i)
when the true genetic model is recessive,(15)(Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′|H1~N6(μR,σ2(V1ρρ′VJ)),whereμR=(μ1RR,μ1RA,μ1RD,μJRR,μJRA,μJRD)′ with
(16)μ1RR=-(n10+n11)n12n1β1,μ1RA=-n12(2n10+n11)β1n1[n10(n11+4n12)+n11n12],μ1RD=-n10n12β1n1(n11+n12),μJRR=-(N0+N1)N2nβ1,μJRA=-N2(2N0+N1)β1n[N0(N1+4N2)+N1N2],μJRD=-N0N2β1n(N1+N2),(ii)
when the true genetic model is additive,(17)(Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′|H1~N6(μA,σ2(V1ρρ′VJ)),whereμA=(μ1AR,μ1AA,μ1AD,μJAR,μJAA,μJAD)′ with
(18)μ1AR=-n12(n11+2n10)β1n1(n10+n11),μ1AA=-n10(n11+4n12)+n11n12n1β1,μ1AD=-n10(n11+2n12)β1n1(n11+n12),μJAR=-N2(N1+2N0)β1n(N0+N1),μJAA=-N0(N1+4N2)+N1N2nβ1,μJAD=-N0(N1+2N2)β1n(N1+N2),(iii)
when the true genetic model is dominant,(19)(Z1R,Z1A,Z1D,ZJR,ZJA,ZJD)′|H1~N6(μD,σ2(V1ρρ′VJ)),whereμD=(μ1DR,μ1DA,μ1DD,μJDR,μJDA,μJDD)′ with
(20)μ1DR=-n12n10β1n1(n10+n11),μ1DA=-n10(n11+2n12)β1n1[n10(n11+4n12)+n11n12],μ1DD=-n10(n11+n12)n1β1,μJDR=-N2N0β1n(N0+N1),μJDA=-N0(N1+2N2)β1n[N0(N1+4N2)+N1N2],μJDD=-N0(N1+N2)nβ1.We develop a method for simplifying the calculations ofPrH0(F1MAX≤u1) and PrH0(FJMAX≤uJ) and PrH0(F1MAX≤u1,FJMAX≤uJ). The details are included in Appendix B, and the calculations of PrH1(F1MAX≤u1) and PrH1(FJMAX≤uJ) and PrH1(F1MAX≤u1,FJMAX≤uJ) are essentially similar.


## 3. Results

### 3.1. Power Comparison

We conduct simulation studies to evaluate the performance of the proposed method under three commonly used genetic models (recessive, additive, and dominant models). We mainly compare the power of two approaches; one is the proposed method in this paper, and the other is the joint analysis based on theF-test statistics F1A and FJA. For convenience, we refer to the proposed method as MAXFJ and AFJ for the other one. We choose the sample size n=2000, and m=5×105. The proportion of subjects genotyped in Stage 1 has three levels π=0.3,0.4,0.5. We set the genome-wide significance level as α=0.05 and that the significance level per SNP as α/m=1×10-7. In Stage 1, the P-value threshold for SNPs selected for followup is set to be 1×10-4 and 2×10-4. We assume that the Hardy-Weinberg equilibrium holds in the general sample population, and then there are on average n×(1-MAF)2, 2n×MAF×(1-MAF), and n×MAF2 individuals with genotype gg, Gg, and GG, respectively, where the minor allele frequency is set to be 0.15, 0.30 and 0.45. To make the power comparison more distinctly, we specify different genetic effect parameters β1 under three genetic models as follows: β1=0.5 for the recessive model, β1=0.3 for the additive model, and β1=0.4 for the dominant model.The power results are displayed in Tables1 and 2 for γ=1×10-4 and γ=2×10-4, respectively. They indicate that MAXFJ is more efficiency robust than AFJ across various inheritance models. As expected, AFJ is more powerful than MAXFJ under the additive model. However, MAFJ performs much more powerful than AFJ when the true genetic model is recessive. For instance, in Table 2, with π=0.4 and MAF = 0.3, the powers of AFJ and MAXFJ are 0.101 and 0.529, respectively. In summary, MAXFJ is substantially more powerful than AFJ in two-staged GWAS of quantitative traits, when the model for AFJ is misspecified.Table 1
Power comparison (n=2000,γ=1×10-4,α=0.05,  and
m=5×105).
π
MAF
REC
ADD
DOM
AFJ
MAXFJ
AFJ
MAXFJ
AFJ
MAXFJ
0.30
0.15
7.5
e
-
5
0.005
0.426
0.365
0.610
0.618
0.30
0.052
0.285
0.811
0.759
0.698
0.784
0.45
0.487
0.785
0.893
0.854
0.449
0.647
0.40
0.15
1.1
e
-
4
0.009
0.651
0.589
0.826
0.837
0.30
0.086
0.470
0.945
0.922
0.887
0.938
0.45
0.711
0.938
0.979
0.968
0.677
0.859
0.50
0.15
1.0
e
-
4
0.010
0.802
0.751
0.933
0.941
0.30
0.121
0.639
0.987
0.980
0.965
0.986
0.45
0.856
0.987
0.997
0.995
0.826
0.953Table 2
Power comparison (n=2000,γ=2×10-4,α=0.05,  and m=5×105).
π
MAF
REC
ADD
DOM
AFJ
MAXFJ
AFJ
MAXFJ
AFJ
MAXFJ
0.30
0.15
1.3
e
-
4
0.006
0.489
0.426
0.676
0.681
0.30
0.066
0.340
0.852
0.806
0.754
0.828
0.45
0.556
0.833
0.922
0.891
0.516
0.706
0.40
0.15
1.2
e
-
4
0.011
0.709
0.651
0.866
0.876
0.30
0.101
0.529
0.961
0.943
0.916
0.956
0.45
0.765
0.957
0.987
0.979
0.732
0.892
0.50
0.15
1.7
e
-
4
0.012
0.838
0.793
0.951
0.958
0.30
0.133
0.683
0.992
0.987
0.975
0.991
0.45
0.888
0.992
0.998
0.997
0.860
0.967


### 3.2. An Illustration Example: Rheumatoid Arthritis

Rheumatoid arthritis (RA) is an autoimmune disease (resulting in a chronically systemic inflammatory disorder) which mainly attacks synovial joints. About 1% of the common adult population worldwide is affected by RA [16]. It has been pointed out that the genetic variants might play a major role in RA susceptibility [17]. Genetic Analysis Workshop 16 (GAW16) based on the North American Rheumatoid Arthritis Consortium (NARAC) is a GWAS testing association with RA using about 5×105 SNPs [18–20]. It included 868 individuals who were RA positive (cases) and also had continuous trait anticyclic citrullinated peptide (anti-CCP) measures and 1194 controls sampled from the New York Cancer Project (NYCP) without RA which had no anti-CCP measures. Huizinga et al. [21] pointed out that a greater anti-CCP would be linked to better prediction of increased risk developing RA. Chen et al. [22] showed that SNP rs2476601 located in PTPN22 had the most significant association with RA. Here, we only focus on SNP rs2476601 and apply two joint analysis methods (AFJ and MAXFJ) to evaluate its statistical significance. The minimum of anti-CCP among 868 cases was affected to each control, and a log transformation of anti-CCP was applied in the analysis. Then, we considered π=0.3,0.4,0.5 three simulation circumstances. For π=0.3, thirty percent of individuals were randomly sampled from all cases and controls and were used as the data from Stage 1, and the rest of individuals were treated as the data of Stage 2. The P-values of AFJ and MAXFJ were calculated, respectively. We repeated the above procedure 1,000 times and saved the corresponding P-values. A base-10 logarithm transformation and an opposite transformation were successively applied to these P-values, and the histogram and density of these transformed data were obtained (Figure 1). Similarly, we conducted the simulation and calculation for π=0.4 and 0.5, and the corresponding histogram and density were presented in Figures 2 and 3. Examination of Figures 1–3 showed that the P-values of MAXFJ are more stable than those of AFJ and the estimated density curves of MAXFJ are more closer to the symmetrical normal distribution while the estimated density curves of AFJ are rather skewed, which indicated that MAXFJ possesses more robust performance when the real genetic models are unknown.The histogram and density of-log10P when π=0.3 (the left subgraph corresponds to MAXFJ while the right one for AFJ).
(a)
(b)The histogram and density of-log10P when π=0.4 (the left subgraph corresponds to MAXFJ while the right one for AFJ).
(a)
(b)The histogram and density of-log10P when π=0.5 (the left subgraph corresponds to MAXFJ while the right one for AFJ).
(a)
(b)


## 3.1. Power Comparison

We conduct simulation studies to evaluate the performance of the proposed method under three commonly used genetic models (recessive, additive, and dominant models). We mainly compare the power of two approaches; one is the proposed method in this paper, and the other is the joint analysis based on theF-test statistics F1A and FJA. For convenience, we refer to the proposed method as MAXFJ and AFJ for the other one. We choose the sample size n=2000, and m=5×105. The proportion of subjects genotyped in Stage 1 has three levels π=0.3,0.4,0.5. We set the genome-wide significance level as α=0.05 and that the significance level per SNP as α/m=1×10-7. In Stage 1, the P-value threshold for SNPs selected for followup is set to be 1×10-4 and 2×10-4. We assume that the Hardy-Weinberg equilibrium holds in the general sample population, and then there are on average n×(1-MAF)2, 2n×MAF×(1-MAF), and n×MAF2 individuals with genotype gg, Gg, and GG, respectively, where the minor allele frequency is set to be 0.15, 0.30 and 0.45. To make the power comparison more distinctly, we specify different genetic effect parameters β1 under three genetic models as follows: β1=0.5 for the recessive model, β1=0.3 for the additive model, and β1=0.4 for the dominant model.The power results are displayed in Tables1 and 2 for γ=1×10-4 and γ=2×10-4, respectively. They indicate that MAXFJ is more efficiency robust than AFJ across various inheritance models. As expected, AFJ is more powerful than MAXFJ under the additive model. However, MAFJ performs much more powerful than AFJ when the true genetic model is recessive. For instance, in Table 2, with π=0.4 and MAF = 0.3, the powers of AFJ and MAXFJ are 0.101 and 0.529, respectively. In summary, MAXFJ is substantially more powerful than AFJ in two-staged GWAS of quantitative traits, when the model for AFJ is misspecified.Table 1
Power comparison (n=2000,γ=1×10-4,α=0.05,  and
m=5×105).
π
MAF
REC
ADD
DOM
AFJ
MAXFJ
AFJ
MAXFJ
AFJ
MAXFJ
0.30
0.15
7.5
e
-
5
0.005
0.426
0.365
0.610
0.618
0.30
0.052
0.285
0.811
0.759
0.698
0.784
0.45
0.487
0.785
0.893
0.854
0.449
0.647
0.40
0.15
1.1
e
-
4
0.009
0.651
0.589
0.826
0.837
0.30
0.086
0.470
0.945
0.922
0.887
0.938
0.45
0.711
0.938
0.979
0.968
0.677
0.859
0.50
0.15
1.0
e
-
4
0.010
0.802
0.751
0.933
0.941
0.30
0.121
0.639
0.987
0.980
0.965
0.986
0.45
0.856
0.987
0.997
0.995
0.826
0.953Table 2
Power comparison (n=2000,γ=2×10-4,α=0.05,  and m=5×105).
π
MAF
REC
ADD
DOM
AFJ
MAXFJ
AFJ
MAXFJ
AFJ
MAXFJ
0.30
0.15
1.3
e
-
4
0.006
0.489
0.426
0.676
0.681
0.30
0.066
0.340
0.852
0.806
0.754
0.828
0.45
0.556
0.833
0.922
0.891
0.516
0.706
0.40
0.15
1.2
e
-
4
0.011
0.709
0.651
0.866
0.876
0.30
0.101
0.529
0.961
0.943
0.916
0.956
0.45
0.765
0.957
0.987
0.979
0.732
0.892
0.50
0.15
1.7
e
-
4
0.012
0.838
0.793
0.951
0.958
0.30
0.133
0.683
0.992
0.987
0.975
0.991
0.45
0.888
0.992
0.998
0.997
0.860
0.967


## 3.2. An Illustration Example: Rheumatoid Arthritis

Rheumatoid arthritis (RA) is an autoimmune disease (resulting in a chronically systemic inflammatory disorder) which mainly attacks synovial joints. About 1% of the common adult population worldwide is affected by RA [16]. It has been pointed out that the genetic variants might play a major role in RA susceptibility [17]. Genetic Analysis Workshop 16 (GAW16) based on the North American Rheumatoid Arthritis Consortium (NARAC) is a GWAS testing association with RA using about 5×105 SNPs [18–20]. It included 868 individuals who were RA positive (cases) and also had continuous trait anticyclic citrullinated peptide (anti-CCP) measures and 1194 controls sampled from the New York Cancer Project (NYCP) without RA which had no anti-CCP measures. Huizinga et al. [21] pointed out that a greater anti-CCP would be linked to better prediction of increased risk developing RA. Chen et al. [22] showed that SNP rs2476601 located in PTPN22 had the most significant association with RA. Here, we only focus on SNP rs2476601 and apply two joint analysis methods (AFJ and MAXFJ) to evaluate its statistical significance. The minimum of anti-CCP among 868 cases was affected to each control, and a log transformation of anti-CCP was applied in the analysis. Then, we considered π=0.3,0.4,0.5 three simulation circumstances. For π=0.3, thirty percent of individuals were randomly sampled from all cases and controls and were used as the data from Stage 1, and the rest of individuals were treated as the data of Stage 2. The P-values of AFJ and MAXFJ were calculated, respectively. We repeated the above procedure 1,000 times and saved the corresponding P-values. A base-10 logarithm transformation and an opposite transformation were successively applied to these P-values, and the histogram and density of these transformed data were obtained (Figure 1). Similarly, we conducted the simulation and calculation for π=0.4 and 0.5, and the corresponding histogram and density were presented in Figures 2 and 3. Examination of Figures 1–3 showed that the P-values of MAXFJ are more stable than those of AFJ and the estimated density curves of MAXFJ are more closer to the symmetrical normal distribution while the estimated density curves of AFJ are rather skewed, which indicated that MAXFJ possesses more robust performance when the real genetic models are unknown.The histogram and density of-log10P when π=0.3 (the left subgraph corresponds to MAXFJ while the right one for AFJ).
(a)
(b)The histogram and density of-log10P when π=0.4 (the left subgraph corresponds to MAXFJ while the right one for AFJ).
(a)
(b)The histogram and density of-log10P when π=0.5 (the left subgraph corresponds to MAXFJ while the right one for AFJ).
(a)
(b)


## 4. Discussion

We have developed a feasible two-stage design and the corresponding robust joint analysis approach for quantitative trait GWASs. The method is based on theF-statistics over three different genetic models. The denominator of the used F-statistic, which is constructed without assuming any genetic model, is different from the commonly used one. This adoption can reduce the computation intensity. Taking advantage of an ingenious design matrix, we successfully construct the common denominator of three F-test statistics for the joint analysis with combined raw data from both stages. The statistical significance (P-value) for the proposed joint analysis method can be calculated with the derived analytic expressions on the basis of the asymptotic distributions, which greatly reduce the complexity and computational intensity compared with the resampling-type permutation and bootstrap procedures. Our numerical results demonstrate that this novel approach has the greater efficiency robustness for genetic model uncertainty than the F-statistic-based joint analysis which assumes the additive genetic model.In this work, we did not investigate the power of joint analysis based on other existing robust association methods for quantitative traits such as So and Sham’s method. We find that it is very difficult to extend So and Sham’s method (score test-based MAX3) to two-staged GWASs with quantitative outcomes, since it is almost impossible to derive the joint distribution of score tests from two stages.For simplicity, here we do not take into account the effects of covariates in the considered two-stage design. However, in real application, the proposed method can be easily applied to the situation including one or more covariates as shown by the original MAXF by Li et al. [10]. It is important to stress that we combine the raw data from two stages to construct the joint statistic, unlike the joint analysis for binary traits using the weighted sum of two statistics in Stages 1 and 2 [12]. Furthermore, one basic assumption in this paper is that the effect sizes of genetic variants between two stages are identical (i.e., no heterogeneity exists), which is the natural and reasonable precondition for the data fusion strategy. In addition, the population-based genetic association studies may be affected by the population stratification, and this needs future research to examine it.



---

*Source: 843563-2013-10-31.xml*