On the Analysis of Genome-Wide Association Studies in Family-Based Designs: A Universal, Robust Analysis Approach and an Application to Four Genome-Wide Association Studies

For genome-wide association studies in family-based designs, we propose a new, universally applicable approach. The new test statistic exploits all available information about the association, while, by virtue of its design, it maintains the same robustness against population admixture as traditional family-based approaches that are based exclusively on the within-family information. The approach is suitable for the analysis of almost any trait type, e.g. binary, continuous, time-to-onset, multivariate, etc., and combinations of those. We use simulation studies to verify all theoretically derived properties of the approach, estimate its power, and compare it with other standard approaches. We illustrate the practical implications of the new analysis method by an application to a lung-function phenotype, forced expiratory volume in one second (FEV1) in 4 genome-wide association studies.


Introduction
During the analysis phase of genome-wide association studies, one is confronted with numerous statistical challenges. One of them is the decision about the ''right'' balance between maximization of the statistical power and, at the same time, robustness against confounding. In family-based designs, the possible range of analysis options spans from a traditional family-based association analysis [1][2][3][4], e.g. TDT, PDT, FBAT, to the application of population-based analysis methods that have been adapted to family-data [1][2][3]. While, by definition, the first group of approaches is completely immune to population admixture and model misspecification of the phenotype, and can be applied to any phenotype that is permissible in the family-based association testing framework (FBAT [4][5][6]), the second category of approaches maximizes the statistical power by a population-based analysis. The phenotypes are modeled as a function of the genotype, and population-based methods such as genomic control [7,8], STRUCTURE [9] and EIGENSTRAT [10], are applied to account for the effects of population admixture and stratification. Hybrid-approaches that combine elements of both populationbased and family-based analysis methods, e.g. VanSteen algorithm [11] and Ionita weighting-schemes [12,13] have been suggested to bridge between the 2 types of analysis strategies. Contrary to the other methods that combine family data and unrelated samples [14][15][16][17], such hybrid testing strategies maintain the 2 key features of the family-based association tests: The robustness against confounding due to population admixture and heterogeneity, and the analysis flexibility of the approach with respect to the choice of the target phenotype. Such 2-stage testing strategies utilize the information about the association at a population-level, the between-family component, to prioritize SNPs for the second step of the approach in which they are tested formally for association with a family-based test. The hybrid approaches can achieve power levels that are similar to approaches in which standard population-based methods are applied to family-data, but the optimal combination of the 2 sources of information (the betweenfamily component and the within-family component) is not straightforward in the hybrid approaches.
In this communication, we propose a new family-based association test for genome-wide association studies that combines all sources of information about association, the between and the within-family information, into one single test statistic. The new test is robust against population-admixture even though both components, the between and the within-family components, are used to assess the evidence for association. The approach is applicable to all phenotypes or combinations of phenotypes that can be handled in the FBAT-approach, e.g. binary, continuous, time-to-onset, multivariate, etc [4][5][6]18]. While the correct model specification for the phenotypes will increase the power of the proposed test statistic, misspecification of the phenotypic model does not affect the validity of the approach. Using extensive simulation studies, we verify the theoretically derived properties of the test statistic, assess its power and compare it with other standard approaches. An application to the Framing heart study (FHS) illustrates the value of the approach in practice. A new genetic locus for the lung-function phenotype, FEV1 (forced expiratory volume in the first second) is discovered and replicated in 3 independent, genome-wide association studies.

Methods
We assume that in a family-based association study, n family members have been genotyped at m loci with a genome-wide SNPchip. For each marker locus, a family-based association test is constructed based on the offspring phenotype and the withinfamily information. The within-family information is defined as the difference between the observed, genetic marker score and the expected, genetic marker score, which is computed conditional upon both the parental genotypes/sufficient statistic [19] under the assumption of Mendelian transmissions. We denote the family-based association test for the ith marker locus by FBAT i . Such an FBAT statistic can be the standard TDT, an FBAT for quantitative/qualitative traits, FBAT-GEE for multivariate traits, etc [4,6,18,20,21]. Similarly, for the ith marker, the betweenfamily information can be used to assess the evidence for association at a population-based level by computing a VanSteen-type [11] ''screening statistic'' T i . The screening statistic is computed based on the data for offspring phenotype and the parental genotypes/sufficient statistic. The statistic T i can be a Wald test for the genetic effect size that is estimated based on the conditional mean model [22], or the estimated power of the family-based test FBAT i [23], either of which is feasible. However, while the FBAT i statistic is robust against population stratification, the screening statistic T i is susceptible to confounding. For this reason, the VanSteen-type testing strategies have restrictively used the between-family information as weights for p-values of the FBAT-statistic, but never as a component in the test statistic itself.

Construction of an overall family-based association test including the population-based and family-based components
In order to construct a family-based association test that incorporates both the within and the between-family information, the Z-statistics that correspond to the p-values of FBAT i and T i are computed. The statistic Z a * is the a quantile of standard normal distribution. pFBAT i and pT i are the p-value of the FBAT-statistics and one sided p-value of the screening statistic where the direction of the one sided screening statistic is defined by the directionality of FBAT i . Based on the statistical independence of FBAT i and T i [11] under the null-hypothesis, we can obtain an overall familybased association test statistic Z i by combining both Z-statistics in a weighted sum, where the parameters w FBAT and w T are standardized weights so that the overall family-based association test Z i has a normal distribution with mean 0 and variance 1, i.e. w FBAT 2 +w T 2 = 1. In the literature, this approach of combining two test statistics is known as the Liptak-method [24]. However, the Liptak-method varies here from its standard application in that the 2 test statistics have to be combined so that confounding in the screening statistic T i cannot affect the validity of the overall family-based association test statistic Z i . In the context of a genome-wide association study (GWAS), we are able to achieve this goal by using rank-based pvalues for the screening statistic T i instead of their asymptotic p-values.
The ''screening statistics'' T i are sorted based on their evidence for association so that T (m) denotes the screening statistic with the least amount of evidence for association and T (1) the screening statistic with the largest amount of evidence for association. The rank-based p-value, (i -0.5)/m, is then assigned to the ordered screening test statistics T (i) . If there is a tie, then the average of the ranks will be used for the computation of the rank-based p-value for the ith marker. Since the null-hypothesis will be true for the vast majority of the SNPs in a GWAS, the rank-based p-values provide an alternative way to assess the significance of the population-based screening statistic T i . The overall association test Z i is then computed based on the Z-score for the asymptotic pvalue of the FBAT-statistic and the Z-score for the ranked-based pvalue of the screening statistic T i . In Text S1 we show that the overall association test Z i maintains the global significance level a, under any situations including population admixture and

Author Summary
In genome-wide association studies, the multiple testing problem and confounding due to population stratification have been intractable issues. Family-based designs have considered only the transmission of genotypes from founder to nonfounder to prevent sensitivity to the population stratification, which leads to the loss of information. Here we propose a novel analysis approach that combines mutually independent FBAT and screening statistics in a robust way. The proposed method is more powerful than any other, while it preserves the complete robustness of family-based association tests, which only achieves much smaller power level. Furthermore, the proposed method is virtually as powerful as populationbased approaches/designs, even in the absence of population stratification. By nature of the proposed method, it is always robust as long as FBAT is valid, and the proposed method achieves the optimal efficiency if our linear model for screening test reasonably explains the observed data in terms of covariance structure and population admixture. We illustrate the practical relevance of the approach by an application in 4 genome-wide association studies.
stratification. This can be understood intuitively as well. The smallest rank-based p-value is 0.5/m. Using the Bonferronicorrection to adjust for multiple testing, the individual, adjusted significance level is a/m which will always be smaller than the smallest rank-based p-value, 0.5/m, unless the pre-specified global significance level a is great than 0.5. This implies that the overall family-based association test can never achieve genome-wide significance just based on the rank-based p-values alone. The FBAT-statistic has to contribute evidence for the association as well in order for the overall family-based association test to reach genome-wide significance. Finally, we have to address the specification of the weights w FBAT and w T in the overall familybased association test statistic Z i . While any combination of weights w FBAT and w T will provide a valid test statistic Z i , the most powerful overall statistic Z i is approximately achieved when the ratio of the weights is equal to the ratio of the standardized effect sizes, the expected effect size of the regression coefficient divided by its (estimated) standard deviation. For quantitative traits in unascertained samples, one can show that optimal power levels are achieved for equal weights, i.e. w FBAT = w T . In general, the equal weighting scheme seems to provide good power levels for any disease mode of inheritance and for different trait types, e.g. binary traits, time-to-onset, etc. The theoretical derivation of optimal weighting schemes for such scenarios is ongoing research and will be published subsequently.
Furthermore, it is important to note that, instead of the Liptakmethod, Fisher's method for combining p-values could have been used as well to construct an overall family-based association test which would have the same robustness properties as the overalltest based on the Liptak-method. However, simulation studies (data not shown) suggest that the highest power levels are consistently achieved with the Liptak method. We therefore omit the approach based on Fisher's method here.

Type I error for 500K GWAS
In the first part of the simulation study, the type-1 error of the proposed family-based association test denoted as LIP was assessed in the absence and in the presence of population admixture, and we use the Wald test based on the conditional mean model [22] with between-family component for pT i in our all simulations. For various scenarios, we verified that the proposed overall familybased association test maintains the a-level.
For simplicity, we assume in the simulation studies that the random samples are given, i.e. no ascertainment, and that the parental genotypes are known. Assuming Hardy-Weinberg equilibrium, the parental genotypes are generated by drawing from Bernoulli distributions defined by the allele frequencies. The offspring genotypes are obtained by simulated Mendelian transmissions from the parents to the offspring. For the jth trio, the offspring phenotype Y j is simulated from a Normal distribution with mean aX j and variance 1, i.e. N(aX j , 1), where the parameter a represents the genetic effect size and the variable X j denotes the offspring genotype. Under the null-hypothesis of no association, the genetic effect size parameter a will be set to 0.
For scenarios in which population admixture is present, we assume that the admixture is created by the presence of 2 subpopulations whose phenotypic means differ by 0.2. The allele frequencies for each marker in the two subpopulations are generated by the Balding-Nichols model [25]. That is, for each marker, the allele frequency in an ancestral population is generated from a uniform distribution between 0.1 and 0.9, U(0.1, 0.9). Then, the marker allele frequencies for the two subpopulations are independently sampled from the beta distri- for the whole markers in each replicate of the simulated GWAS. A survey reported F ST estimates with a median of 0.008 and 90th percentile of 0.028 among Europeans, and the corresponding values are 0.027 and 0.14 among Africans, and 0.043 and 0.12 among Asians [26]. The value for Wright's F ST was assumed to be 0.05, 0.1, 0.2, or 0.3. Each trio was assigned to the one of the 2 subpopulations with 50% probability.
In the absence and presence of the population stratification (F ST = 0.05, 0.1, 0.2, and 0.3), Table 1 shows the empirical type-1 error rates of the overall association test statistic Z i for a GWAS with 500,000 SNPs. The estimates for the empirical significance levels in Table 1 are based on 2,000 replicates. The empirical genome-wide significance level is calculated as the proportion of replicates for which the minimum p-values among the 500,000 markers is less than 0.05/500,000. We consider the proposed equal weights for w FBAT and w T , for genome-wide significance level 0.05 and Table 1 shows that the type-1 error rate is preserved well. For different significance levels, we calculate in Table 2 the empirical proportions of SNPs for which the overall family-based association test Z i is significant at the a-levels of 0.05, 0.01, 10 23 , 10 24 and 10 25 . The simulation studies are conducted in the absence and in the presence of population admixture. Table 2 does not provide any evidence for a departure of the empirical significance levels from the theoretical levels, both in the absence and presence of population substructure. These results confirm our theoretical conclusions that Z i is robust against population stratification and maintains correct type-1 error.
In the next set of simulation studies, we assess the effects of the local population stratification on the overall family-based association test. We generate local population stratification under the  The number of trios, N trio , is assumed to be 1,000 and the empirical proportions of SNPs whose p-values for Z i are less than c in each replicate for 500K GWAS are averaged over 2,000 replicates. doi:10.1371/journal.pgen.1000741.t002 following assumptions: there are two subpopulations, G 1 and G 2 which distinguish themselves from each other in 2 marker regions. We assume that a subject can be from all possible 4 combinations at the 2 particular regions, e.g. (G 1 , G 1 ), (G 1 , G 2 ), (G 2 , G 1 ) and (G 2 , G 2 ). Both regions consist of 10K SNPs and 90K SNPs respectively and if subjects are from the same subpopulation in each genetic region, their assumed allele frequencies of the markers in the corresponding region are equal. For example, the allele frequencies of each marker in the marker region 1 are the same for samples in (G 1 , G 1 ) and (G 1 , G 2 ), but they are different for (G 1 , G 1 ) and (G 2 , G 2 ). In the simulation study, we generate the parental genotypes based on these allele frequency assumptions and obtain the offspring genotypes based on simulated Mendelian transmissions. Using the Balding-Nichols model we considered F ST 's of 0.001, 0.005, 0.01 and 0.05 in the simulation studies. The offspring's phenotype was generated under the null hypothesis, but we assumed that each sub-population strata had a different phenotypic mean: 0 for (G 1 , G 1 ), 0.2 for (G 1 , G 2 ), 0.4 for (G 2 , G 1 ) and 0.6 for (G 2 , G 2 ). Each replicate consists of 2,000 trios with equal number of trios for all 4 possible combinations. The data was analyzed with the proposed overall family-based association test and with standard linear regression after adjusting population admixture with EIGENSTRAT [10]. For EIGENSTRAT, we applied the principal component analysis to the mean of the paternal and maternal genotypes at each locus because parents of each offspring are from the same subpopulation, and then the residuals obtained from regressing offspring genotypes and phenotypes with eigenvectors respectively are used to calculate the generalized Armitage trend test [27]. Table 3 provides the empirical type-1 error for both analysis approaches based on 2,000 replicates. While EIGENSTRAT exhibits an inflated type-1 error, the proposed overall family test maintains the theoretical significance level.
Empirical power with simulation for 500K GWA for quantitative trait For the analysis of quantitative traits, Table 4 provides the empirical power for 500K GWAS from 2000 replicates when there is no population stratification. Under the assumption of an additive disease model for a quantitative trait, the genetic effect, a, is given as a function of the heritability, h 2 , the minor allele frequency p Dı and the phenotypic variance, s 2 , by: a = sh/ [2p(12p)(12h 2 ) ] 0.5 . In the simulation study, we assume heritabilities of h 2 = 0.001, 0.005, 0.01 and 0.015 for 2,000, 2,500 and 3,000 trios. The allele frequency of the disease locus, p Dı , is 0.3 and the phenotypic variance is 1. We compare the achieved power levels of the proposed overall family-based association test, Z i , with the weighting approach by Ionita-Laza et al [12], the original VanSteen approach [11], the QTDT approach [28] and population-based analysis, i.e. using linear regression of the phenotype Y on the genotype X. Bonferroni correction is used to adjust for multiple testing in the population-based analysis, FBAT, QTDT and the proposed method. The results in Table 4 suggest that the proposed association test achieves power levels that represent a major improvement over the existing methods for family-based association tests (VanSteen [11] or Ionita-Laza [12]). Our approach reaches the same power levels as the populationbased analysis. For the power comparisons that are shown in Figure 1, Figure 2, and Figure 3, the number of trios is assumed to be 1,000 in 500K GWAS and the empirical powers are calculated based on 10,000 replicates at an a-level of 0.001 for the all genetic  The number of trios, N trio , is assumed to be 1,000. The empirical proportions of SNPs whose p-values for Z i are less than c in each replicate for 500K GWAS are averaged over 2000 replicates when there is local population stratification. LIP stands for the proposed method using Liptak method to combine pFBAT i and pT i . doi:10.1371/journal.pgen.1000741.t003 models. The results confirm that the Liptak's method combining T i and FBAT i has similar power to the population-based method, and the choice of equal weights performs well. The simulation results in Table 4 also suggest that QTDT [28] approach achieves similar power levels as the standard FBAT approach, which is consistent with previously reported findings in the literature [29]. However, both standard FBAT and QTDT are still much less powerful than the proposed overall family-based association test. Table 5 shows the empirical power for a GWAS with 100,000 SNPs in the presence of population stratification. For the parameters of this simulation study, we assume F ST = 0.001, 0.005, 0.01, and 0.05, and the additive mode of inheritance at the disease locus with values for the heritability of h 2 = 0.005, 0.01 and 0.015. The disease allele frequency p Dı in the ancestral population is assumed to be 0.3. The phenotypic data is simulated so that their phenotypic means for two subpopulations are 0 and 0.2 respectively. Each individual/trio is assigned to either subpopulation with probability 0.5. The parental genotypes are used to estimate the ancestry for EIGENSTRAT as before. Various methods have been suggested to adjust the population stratification in a population-based designs and we compare the proposed methods with the EIGENSTRAT approach [10]. In order to maximize the power of the proposed method, we apply the EIGENSTRAT approach to the population-based component T i of our approach, i.e. principal component analysis based on the parental genotypes and the offspring's phenotype is integrated into the generalized Armitage test for T i [27]. To keep the power comparisons unbiased, the population-based components of the approaches by VanSteen and Ionita-Laza are also adjusted for population admixture, using the EIGENSTRAT approach. The results in Table 5 show that the proposed test statistic Z i is considerably more powerful than population-based analysis adjusted with EIGENSTRAT. QTDT is slightly more powerful than FBAT, but it is much less powerful than LIP as is in Table 4. This suggests that EIGENSTRAT should be applied only to between-family component in family-based association studies.
Our unpublished work showed that the proposed approach can be less powerful than the combination of population-based analysis and EIGENSTRAT if pT i is calculated from the conditional mean model [11,22] without adjusting population stratification.

Applications to a genome-wide association in the Framingham Heart study
For the assessment of the severity of pulmonary diseases, the lung volume of air that a subject can blow out within one second after taking a deep breath is an important endo-phenotype. It is referred to as the forced expiratory volume in one second (FEV1). FEV1 is an important measure for lung function and we apply the proposed method to a family-based GWAS of FEV1. The proposed method is applied to 550K GWAS Framingham Heart Study (FHS) data set for FEV1, and then we confirm whether the selected SNPs are replicated in the British 1958 Birth Cohort (BBC), another population sample, as well as two samples of asthmatics in the the Childhood Asthma management program (CAMP) [30] and an Afro-Caribbean group of families from Barbados (ACG) [31]. In FHS, 9,274 subjects were genotyped and 10,816 subjects of those had at least one FEV1 measurement. Of the 8637 participants with genotyping and FEV1 measures, only those with a call rate of 97% or higher were included. We adjusted the covariates, age, sex and the quadratic term of height that are known to be associated with FEV1. For within-family components, the FBAT statistic for quantitative trait was applied. Markers were excluded from the analysis if the number of informative families was less than 20, or the minor allele frequency was less than 0.05. In total, 306,264 SNPs were used for analysis and, based on the number of SNPs, rankbased empirical p-values, pT i , and the genome-wide significance level was obtained with Bonferroni correction. When we let n and n inf be the total number of individuals and the number of informative trios respectively, n inf : (2n2n inf ) are used for the weights of Z i because some of parental phenotypes are available. Table 6 shows the p-values for the top 10 SNPs from the proposed method. In our analysis, the genome-wide significance level at 0.05 is 1.636610 27 and our results show that only the first ranked SNP, rs805294, is significant at the genome-wide level 0.2 with Bonferroni correction. For rs805294, we also checked the significance in other data sets, BBC, CAMP [30] and ACG [31]. In CAMP, 1215 subjects in 422 families were genotyped and there are 488 informative trios for rs809254 and in ACG, there were only 33 informative trios (Table 7). In the BBC, 1372 unrelated subjects were genotyped with the Affymetrix chip and 1323 unrelated subjects genotyped with the Illumina chip. In CAMP and ACG, age, sex and the quadratic terms of heights were adjusted and in the BBC, age, sex, height, recent chest infection and nurse were adjusted. Table 7 also shows that rs805294 is significant and their directions are same for the considered studies except for the ACG sample. In particular, in the ACG study, the MAF of the SNP is different from other studies, which indicates a different local LD structure; The ACG sample is from an Afro-Caribbean population, contrary to the other studies which only include Caucasian study subjects. In addition, the ACG sample lacks statistical power for this particular SNP, i.e. there are only 33 informative trios in this sample. Thus, the inconsistent finding in the ACG study could be attributable to genetic heterogeneity, i.e. different local LD structure/flip-flop phenomena [32], or insufficient statistical power. For meta analysis, the sample sizes are used as weights for Liptak's method and we use

Discussion
Genome-wide association studies have become one of the most important tools for the identification of new disease loci in the human genome. However, even though advances in genotyping technology have enabled a new generation of genetic association studies that provide robust and replicable findings, population stratification/genetic heterogeneity and the multiple testing problems continue to be the major issues in the statistical analysis that have to be resolved in each study. While family-based association tests provide analysis results that are completely robust against confounding due to population-substructures, the analysis approach is not optimal in terms of statistical power. Numerous approaches have been suggested to minimize this disadvantage of family-based association tests but the previous approaches had to compromise either in terms of robustness or in terms of efficiency.
In this communication, we develop an approach that efficiently utilizes all available data, while maintaining complete robustness against confounding due to population substructure. The proposed methods combines the p-values of the family-based tests (the within-component) with the rank-based p-values for populationbased analysis (the between component) to achieve optimal power levels. The use of rank-based p-values for the population-based component is similar in spirit to the genomic control approach. In principle, the genomic control functions as rescaling the variance inflated due to population stratification under the assumption of the constant F ST . Rank-based p-value directly rescales the statistics based on their ranks, which always generates the uniformly distributed p-value and provides validity even for varying F ST due to local population stratification etc.
Although our simulations are limited to independent unascertained samples and quantitative traits, the proposed work can be easily extended to ascertained samples, large pedigree, or different trait types, etc. By replacing the parental genotypes with the sufficient statistics by Rabinowitz&Laird [19], the FBAT-statistic and the screening-statistic can be adopted straight-forwardly to designs with extended pedigrees [23]. Similarly, parental phenotypes can be incorporated into the conditional mean model [23] or its non-parametric extensions [33] as additional outcome variables. The optimal weights can vary between the different scenarios and further theoretical investigation is currently ongoing, but limited initial simulation studies suggest that equal weights, while not always the most powerful choice in such situation, will always result in more powerful analysis than currently used methods.