Volume 8 Supplement 1
Accounting for relatedness in family-based association studies: application to Genetic Analysis Workshop 18 data
© Eu-ahsunthornwattana et al.; licensee BioMed Central Ltd. 2014
Published: 17 June 2014
In the last few years, a bewildering variety of methods/software packages that use linear mixed models to account for sample relatedness on the basis of genome-wide genomic information have been proposed. We compared these approaches as implemented in the programs EMMAX, FaST-LMM, Gemma, and GenABEL (FASTA/GRAMMAR-Gamma) on the Genetic Analysis Workshop 18 data. All methods performed quite similarly and were successful in reducing the genomic control inflation factor to reasonable levels, particularly when the mean values of the observations were used, although more variation was observed when data from each time point were used individually. From a practical point of view, we conclude that it makes little difference to the results which method/software package is used, and the user can make the choice of package on the basis of personal taste or computational speed/convenience.
A number of different methods/software packages have been proposed in the last few years that implement linear mixed-model approaches to account for population structure and relatedness among samples in genome-wide association studies (GWAS), but no detailed comparisons among them have been made before our effort. Indeed, when a new method/package is developed, it is often quite unclear whether or how it differs substantially from those already available. To address this question, we explored the performance of various implementations of such methods in the longitudinal Genetic Analysis Workshop 18 (GAW18) data set.
We analyzed the GAW18 GWAS data  using the real phenotypes and the first set of simulated phenotypes. This analysis was performed without knowledge of the underlying simulating model. The genotype data were cleaned using standard procedures . This resulted in 4 individuals being excluded because of their total lack of genotype data, and another individual being excluded because of outlying ethnicity (Chinese [CHB] or Japanese [JPT]), leaving 954 individuals whose genotype data were used. We removed 43,987 monomorphic or low-frequency (minor allele frequency [MAF] <1%) single-nucleotide polymorphisms (SNPs), 109 SNPs with missing rate above 10% (this criterion took into account the apparently high missing rate in some SNPs likely to be caused by the differences in genotyping technology used in the samples), and 1 SNP that failed Hardy-Weinberg equilibrium testing in the control founder population. A total of 427,952 SNPs were retained for analysis.
We conducted linear regression of the real and simulated systolic blood pressure and simulated diastolic blood pressure at each time point regressed on age, medication, and smoking status. For the real diastolic blood pressure--which, as could be physiologically expected, seemed to have a nonlinear relationship with age--we used a quadratic regression, including age and age squared as predictors. The phenotype data from all individuals were used for these regressions. Residuals from these regressions in subjects who also had genotype data were then used for the genome-wide analyses.
Genome-wide association analyses, adjusting for familial relatedness using genomic data, were performed using a variety of linear mixed model approaches. All approaches attempt to fit the model Y=β+Q+ε, where Y=(y1, ..., y n ) T is a vector of responses on n subjects; X= (x ik ) is the n × K matrix of predictor values for variables to be modeled as fixed effects (including covariates and genotypes at any SNPs currently under test); β=(β1, ... βK) T are regression coefficients (to be estimated) representing the linear effects of the predictors on the response; Q are random effects, Q~N(0,2σg2Φ), and ε are random errors, ε~N(0,σe2I), where σg2 and σe2 are parameters (to be estimated) representing the genetic and environmental components of variance respectively; Φ is the n × n matrix of pairwise kinship coefficients; and I is the n × n identity matrix. The approaches vary with respect to precise details of the calculation of kinship or "relatedness" and with respect to whether an exact method or a fast approximation is used (for more details, see descriptions in references [3–9]). In each case we used a subset of 21,153 SNPs to perform the relatedness calculations, namely SNPs with MAF >0.4, <5% missing data, and "pruned" to be in approximate linkage equilibrium via the PLINK command "-indep 50 5 2". In analyses of other data sets we have found little difference between results when using such a pruned set of SNPs for calculating relatedness and when using the full set of SNPs (data not shown).
The methods considered were: (a) EMMAX , which implements 2 methods for relatedness calculations: one based on identity-by-state (IBS) sharing and one based on the Balding-Nichols method ; (b) FaST-LMM , which also implements 2 methods to adjust for relatedness: one using a standard covariance matrix and one using the realized relationship matrix; (c) the polygenic/mmscore functions in GenABEL , which implement the FASTA method ; (d) the polygenic/grammar functions in GenABEL, which implement the GRAMMAR-Gamma approximation ; and (e) Gemma , which uses an efficient exact method. Simple linear regression without any relatedness adjustment was also performed in FaST-LMM. All analyses were performed using both the residual from each individual observation (modeled without regard to its true longitudinal nature, or longitudinal) and the mean of the residuals for each subject, or mean. Genomic inflation factors (λ) were calculated as proposed by Devlin and Roeder . We also assessed the genomic inflation factors for unadjusted χ2 and Cochran-Armitage trend tests of hypertension status at each time point as calculated using PLINK .
Results and discussion
For the analyses using hypertension status, the unadjusted genomic inflations were between 1.21 and 1.55 for the Cochran-Armitage trend test and between 1.01 and 1.27 for the χ2 test.
Although the results from all packages considered here were similar, the implementations did vary in speed. All packages performed the analysis in reasonable time (less than 1 day) on our system. Precise timings will depend on the computer resources and architecture available, but as a rule of thumb we found FaST-LMM and GRAMMAR-Gamma to be the fastest (taking just a few hours), followed by EMMAX and Gemma, which took 12 to 16 hours, and GenABEL/FASTA, which took 18 to 20 hours.
All methods performed well and results were similar, particularly at the most significant SNPs. We conclude that (at least for nonlongitudinal traits) it makes little difference to the results which method/software package is used, and the user can make the choice of package on the basis of personal taste, speed, or computational convenience. For longitudinal traits (modeled without regard to their longitudinal nature) the slight differences seen between the methods would be an interesting topic for further investigation, but it is beyond the scope of the current article.
This work was supported by the Wellcome Trust (grant reference 087436). JE receives scholarship and funding from Faculty of Medicine, Ramathibodi Hospital, Mahidol University, Bangkok, Thailand. The GAW18 whole genome sequence data were provided by the T2D-GENES Consortium, which is supported by NIH grants U01 DK085524, U01 DK085584, U01 DK085501, U01 DK085526, and U01 DK085545. The other genetic and phenotypic data for GAW18 were provided by the San Antonio Family Heart Study and San Antonio Family Diabetes/Gallbladder Study, which are supported by NIH grants P01 HL045222, R01 DK047482, and R01 DK053889. The Genetic Analysis Workshop is supported by NIH grant R01 GM031575.
This article has been published as part of BMC Proceedings Volume 8 Supplement 1, 2014: Genetic Analysis Workshop 18. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcproc/supplements/8/S1. Publication charges for this supplement were funded by the Texas Biomedical Research Institute.
- Almasy L, Dyer T, Peralta J, Jun G, Fuchsberger C, Almeida M, Kent JW, Fowler S, Duggirala R, Blangero J: Data for Genetic Analysis Workshop 18: human whole genome sequence, blood pressure, and simulated phenotypes in extended pedigrees. BMC Proc. 2014, 8 (suppl 2): S2-PubMed CentralView ArticlePubMedGoogle Scholar
- Anderson CA, Pettersson FH, Clarke GM, Cardon LR, Morris AP, Zondervan KT: Data quality control in genetic case-control association studies. Nat Protoc. 2010, 5: 1564-1573. 10.1038/nprot.2010.116.PubMed CentralView ArticlePubMedGoogle Scholar
- Kang HM, Sul JH, Service SK, Zaitlen NA, Kong SY, Freimer NB, Sabatti C, Eskin E: Variance component model to account for sample structure in genome-wide association studies. Nat Genet. 2010, 42: 348-354. 10.1038/ng.548.PubMed CentralView ArticlePubMedGoogle Scholar
- Balding DJ, Nichols RA: A method for quantifying differentiation between populations at multi-allelic loci and its implications for investigating identity and paternity. Genetica. 1995, 96: 3-12. 10.1007/BF01441146.View ArticlePubMedGoogle Scholar
- Listgarten J, Lippert C, Kadie CM, Davidson RI, Eskin E, Heckerman D: Improved linear mixed models for genome-wide association studies. Nat Methods. 2012, 9: 525-526. 10.1038/nmeth.2037.PubMed CentralView ArticlePubMedGoogle Scholar
- Aulchenko YS, Ripke S, Isaacs A, van Duijn CM: GenABEL: An R library for genome-wide association analysis. Bioinformatics. 2007, 23: 1294-1296. 10.1093/bioinformatics/btm108.View ArticlePubMedGoogle Scholar
- Chen WM, Abecasis GR: Family-based association tests for genomewide association scans. Am J Hum Genet. 2007, 81: 913-926. 10.1086/521580.PubMed CentralView ArticlePubMedGoogle Scholar
- Svishcheva GR, Axenovich TI, Belonogova NM, van Duijn CM, Aulchenko YS: Rapid variance components-based method for whole-genome association analysis. Nat Genet. 2012, 44: 1166-1170. 10.1038/ng.2410.View ArticlePubMedGoogle Scholar
- Zhou X, Stephens M: Genome-wide efficient mixed-model analysis for association studies. Nat Genet. 2012, 44: 821-824. 10.1038/ng.2310.PubMed CentralView ArticlePubMedGoogle Scholar
- Devlin B, Roeder K: Genomic control for association studies. Biometrics. 1999, 55: 997-1004. 10.1111/j.0006-341X.1999.00997.x.View ArticlePubMedGoogle Scholar
- Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, et al: PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007, 81: 559-575. 10.1086/519795.PubMed CentralView ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.