|Home | About | Journals | Submit | Contact Us | Français|
Conceived and designed the experiments: SW CL. Performed the experiments: SW. Analyzed the data: SW JBW RM CJO EKS KB GTO STW. Contributed reagents/materials/analysis tools: SW. Wrote the paper: SW CL.
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.
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 population-based 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.
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 –, e.g. TDT, PDT, FBAT, to the application of population-based analysis methods that have been adapted to family-data –. 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 –), 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 ,, STRUCTURE  and EIGENSTRAT , are applied to account for the effects of population admixture and stratification. Hybrid-approaches that combine elements of both population-based and family-based analysis methods, e.g. VanSteen algorithm  and Ionita weighting-schemes , have been suggested to bridge between the 2 types of analysis strategies. Contrary to the other methods that combine family data and unrelated samples –, 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 between-family 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 –,. 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.
We assume that in a family-based association study, n family members have been genotyped at m loci with a genome-wide SNP-chip. For each marker locus, a family-based association test is constructed based on the offspring phenotype and the within-family 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  under the assumption of Mendelian transmissions. We denote the family-based association test for the ith marker locus by FBATi. Such an FBAT statistic can be the standard TDT, an FBAT for quantitative/qualitative traits, FBAT-GEE for multivariate traits, etc ,,,,. Similarly, for the ith marker, the between-family information can be used to assess the evidence for association at a population-based level by computing a VanSteen-type  “screening statistic” Ti. The screening statistic is computed based on the data for offspring phenotype and the parental genotypes/sufficient statistic. The statistic Ti can be a Wald test for the genetic effect size that is estimated based on the conditional mean model , or the estimated power of the family-based test FBATi , either of which is feasible. However, while the FBATi statistic is robust against population stratification, the screening statistic Ti 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.
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 FBATi and Ti are computed. The statistic Zα* is the α quantile of standard normal distribution. pFBATi and pTi 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 FBATi. Based on the statistical independence of FBATi and Ti  under the null-hypothesis, we can obtain an overall family-based association test statistic Zi by combining both Z-statistics in a weighted sum,
where the parameters wFBAT and wT are standardized weights so that the overall family-based association test Zi has a normal distribution with mean 0 and variance 1, i.e. wFBAT 2+wT 2=1. In the literature, this approach of combining two test statistics is known as the Liptak-method . 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 Ti cannot affect the validity of the overall family-based association test statistic Zi. In the context of a genome-wide association study (GWAS), we are able to achieve this goal by using rank-based p-values for the screening statistic Ti instead of their asymptotic p-values.
The “screening statistics” Ti 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 Ti. The overall association test Zi is then computed based on the Z-score for the asymptotic p-value of the FBAT-statistic and the Z-score for the ranked-based p-value of the screening statistic Ti. In Text S1 we show that the overall association test Zi maintains the global significance level α, under any situations including population admixture and stratification. This can be understood intuitively as well. The smallest rank-based p-value is 0.5/m. Using the Bonferroni-correction to adjust for multiple testing, the individual, adjusted significance level is α/m which will always be smaller than the smallest rank-based p-value, 0.5/m, unless the pre-specified global significance level α 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 wFBAT and wT in the overall family-based association test statistic Zi. While any combination of weights wFBAT and wT will provide a valid test statistic Zi, the most powerful overall statistic Zi 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. wFBAT=wT. 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 Liptak-method, 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 overall-test 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.
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  with between-family component for pTi in our all simulations. For various scenarios, we verified that the proposed overall family-based association test maintains the α-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 Yj is simulated from a Normal distribution with mean aXj and variance 1, i.e. N(aXj, 1), where the parameter a represents the genetic effect size and the variable Xj 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 . 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 distributions (p(1−FST)/FST, (1−p)(1−FST)/FST) for the whole markers in each replicate of the simulated GWAS. A survey reported FST 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 . The value for Wright's FST 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 (FST=0.05, 0.1, 0.2, and 0.3), Table 1 shows the empirical type-1 error rates of the overall association test statistic Zi 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 wFBAT and wT, 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 Zi is significant at the α-levels of 0.05, 0.01, 10−3, 10−4 and 10−5. 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 Zi 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 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 FST'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 . 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 . 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.
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 pD ı and the phenotypic variance, σ 2, by: a=σh/[2p(1−p)(1−h 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, pD ı, is 0.3 and the phenotypic variance is 1. We compare the achieved power levels of the proposed overall family-based association test, Zi, with the weighting approach by Ionita-Laza et al , the original VanSteen approach , the QTDT approach  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  or Ionita-Laza ). Our approach reaches the same power levels as the population-based 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 α-level of 0.001 for the all genetic models. The results confirm that the Liptak's method combining Ti and FBATi 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  approach achieves similar power levels as the standard FBAT approach, which is consistent with previously reported findings in the literature . 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 FST=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 pD ı 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 . In order to maximize the power of the proposed method, we apply the EIGENSTRAT approach to the population-based component Ti 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 Ti . 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 Zi 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 pTi is calculated from the conditional mean model , without adjusting population stratification.
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)  and an Afro-Caribbean group of families from Barbados (ACG) . 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, rank-based empirical p-values, pTi, and the genome-wide significance level was obtained with Bonferroni correction. When we let n and ninf be the total number of individuals and the number of informative trios respectively, ninf: (2n−ninf) are used for the weights of Zi 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.636×10−7 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  and ACG . 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 , or insufficient statistical power. For meta analysis, the sample sizes are used as weights for Liptak's method and we use 131351=FHSBBCCAMPACG as weights because the between-family information is used only for FHS. If the p-value from Illumina gene chip in BBC and the p-values from FHS, CAMP and ACG are combined, then the p-values by Liptak's method using proposed weights and Fisher's method are 1.534×10−8 and 1.081×10−7 respectively, and they become 4.625×10−9 and 3.554×10−8 if the p-values from one-tailed tests are used for BBC, CAMP and ACG with the same direction of FHS. If the p-value from the Affymetrix gene chip in BBC is combined with the other studies, then they are 3.787×10−8 (Liptak's method) and 1.890×10−7 (Fisher's method) for two-tailed tests, and 1.098×10−8 (Liptak's method) and 6.236×10−8 (Fisher's method) for one-tailed tests. As a result we can conclude that rs805294 is significantly associated with FEV1 at a genome-wide scale and the gene, LY6G6C, associated with rs805293 will be investigated in further studies.
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 population-based 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 FST. 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 FST 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 , the FBAT-statistic and the screening-statistic can be adopted straight-forwardly to designs with extended pedigrees . Similarly, parental phenotypes can be incorporated into the conditional mean model  or its non-parametric extensions  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.
The validity of the proposed method.
(0.04 MB DOC)
Framingham Heart Study genotype and phenotype data are publicly available through the NHLBI's SNP Health Association Resource (SHARe) initiative (http://public.nhlbi.nih.gov/GeneticsGenomics/home/share.aspx). We acknowledge the CAMP investigators and research team for collection of CAMP Genetic Ancillary Study data and use of genotype data from the British 1958 Birth Cohort DNA collection. We further acknowledge the families in Barbados for their generous participation in this study. We are grateful to Drs. Raana Naidu, Paul Levett, Malcolm Howitt and Pissamai Maul, Trevor Maul, and Bernadette Gray for their contributions in the field; Dr. Malcolm Howitt and the Polyclinic and A&E Department physicians in Barbados for their efforts and their continued support; as well as Drs. Henry Fraser and Anselm Hennis at the Chronic Disease Research Centre.
The authors have declared that no competing interests exist.
The CAMP Genetics Ancillary Study is supported by U01 HL075419, U01 HL65899, P01 HL083069, R01 HL086601, and T32 HL07427 from the National Heart, Lung, and Blood Institute, National Institutes of Health. CL is supported by the National Institutes of Health grant R01MH081862. Framingham Heart Study genotype and phenotype data are publicly available through the NHLBI's SNP Health Association Resource (SHARe) initiative (http://public.nhlbi.nih.gov/GeneticsGenomics/home/share.aspx). The British 1958 Birth Cohort DNA collection is funded by the Medic Research Council grant G00000934 and the Wellcome Trust grant 068545/Z/02. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.