Adjusting for population stratification and relatedness with sequencing data
© Zhang and Pan; licensee BioMed Central Ltd. 2014
Published: 17 June 2014
To avoid inflated type I error and reduced power in genetic association studies, it is necessary to adjust properly for population stratification and known/unknown subject relatedness. It would be interesting to compare the performance of a principal component-based approach with a linear mixed model. Furthermore, with the availability of genome-wide sequencing data, the question of whether it is preferable to use common variants or rare variants for such an adjustment remains largely unknown. In this paper, we use the Genetic Analysis Workshop 18 data to empirically investigate these issues. We consider both a quantitative trait and a binary trait.
In genetic association studies, population stratification and known or cryptic relatedness are always issues. If not suitably accounted for, these may cause inflated type I errors and reduced power. A popular approach to adjusting for population stratification is to construct principal components (PCs) from some similarity matrix for the samples and to include the PCs as covariates in a regression model . This is referred to as a PC-based approach. However, it is thought that these approaches "do not model family structure or cryptic relatedness" . A more general, and perhaps more powerful, approach is to apply linear mixed models (LMMs) to account for both population stratification and relatedness [3, 4]. EMMAX software  has facilitated the implementation of LMM by using the identity-by-state (IBS) matrix to capture the complex correlation structure in the samples. As these methods are studied intensely in genome-wise association studies (GWAS), a natural question is how they will perform with sequencing data. It is also of interest to investigate whether common variants (CVs) with minor allele frequencies (MAFs) no less than 0.05, or rare variants (RVs) with 0< MAF < 0.01, should be used to infer the samples' genetic similarities.
In this paper, we compare the PC-based approach with LMM to determine which approach can better control the inflation of type I error arising from correlated samples. The association testing is carried out for a quantitative trait and a binary trait in the Genetic Analysis Workshop 18 (GAW18) family-based sequencing data. For a complete comparison, we construct and consider PCs from different similarity matrices: the sample covariance matrix and the IBS matrix. Finally, we discuss the best choice of variants for constructing the similarity matrix, which has been the subject of several recent studies [6–8].
In the PC-based approach, for a given similarity matrix, we obtain its m largest eigenvalues λ j and the corresponding eigenvectors v j for j = 1, ..., m and denote . For a quantitative trait, we use a linear regression model Y = β0 + X m γ + Zζ + gβ + δ , where δ ~ N (0, σ2I). For a binary trait we adopt a logistic model: Logit(E(Y )) = β0 + X m γ + Zζ + gβ.
In the 2 models above, is the vector of the traits for n subjects. is the matrix of covariates, and is the vector of the genotype scores for 1 or more variants to be tested. We denote the method by which PCs are obtained from the IBS matrix as PCA. IBS, and the method by which PCs are obtained from the covariance matrix as PCA.V.
In an LMM, Y = β0 + Zζ + gβ + u + δ, where Y , Z, and g are defined as above, u is the random effect for other polygenic effects and δ is the residual error. It is assumed that δ ~ N(0, σ2I) and , where K is the IBS matrix. We use EMMAX  for parameter estimation and inference.
For hypothesis testing, PCA.V and PCA.IBS adopted the Wald test and EMMAX adopted the F-test.
We used the GAW18 sequencing data containing 959 individuals and 8,348,674 single nucleotide variants (SNVs) across all 11 chromosomes, among which 2,791,923 SNVs were CVs and 3,977,003 were RVs. After pruning by PLINK  using a sliding window of size 50, moving step of 5 and r2 <5%, and filtering out those with missing call >0.05, there were 63,157 CVs left. We randomly selected 10,837 CVs from those to construct the similarity matrix. The IBS matrix was obtained by EMMAX, and the covariance matrix was obtained by the R function cov() with "use = pairwise.complete.obs" to utilize the maximum number of variants.
We used the measurements of systolic blood pressure (SBP) at time point 1, SBP1, and the hypertension diagnosis at time point 1, HTN1. The former is a quantitative trait and the latter is a binary trait. There are 855 samples available. Gender, smoking, and age are the covariates.
Association test with CVs
We carried out single single-nucleotide polymorphism (SNP) analysis on a set of 6228 CVs randomly selected from all the pruned CVs. Based on the findings of previous GWAS that most of the SNPs were not significantly associated with hypertension, we could assume these 6228 CVs were null SNPs. Because some subjects were from the same families and thus correlated, we expected to observe an inflated type I error if we treated the samples as independent. If the PC-based method or LMM was effective in adjustment, the p values should have followed a uniform distribution. This also meant that the proportion of the tests with p value <0.05 should be close to 0.05 and the inflation factor λ close to 1; λ is the inflation factor of p values estimated by the function gcontrol2 in R package gap. It is calculated as the ratio of the medians of the observed and expected statistics, respectively.
Summary statistics of p values for SBP1 by PCA.V, PCA.IBS, and EMMAX. The similarity matrix is based on CVs.
% (pval <0.05)
Summary statistics of p values for HTN1 by PCA.V, PCA.IBS, and EMMAX
% (pval <0.05)
CVs or RVs?
Lastly, we examine which type of variants, CVs or RVs, are more capable of capturing the underlying sample structure. For this purpose, we use PLINK to randomly select 11,103 variants from 1,104,098 pruned RVs to construct the covariance matrix or IBS matrix.
Results of the association tests by PCA
% (p val <0.05)
Originally, the weaker performance of PCA.V based on RVs was thought to be a result of insufficient inclusion of PCs. Following the suggestion of Patterson et al , we use the Tracy-Widom test to test how many PCs are necessary to be considered significant [2, 8]. The test shows that the top 210 PCs of the covariance matrix all have p values smaller than 0.05. However, we fail to obtain reasonable p values with 200 PCs included. This might be because the model could not be fitted, given the small sample size.
In this paper, we address 3 questions: (a) how the PC-based approach and LMM perform in controlling type I error for correlated samples; (b) whether the IBS or covariance matrix should be used to generate PCs; and (c) whether CVs or RVs should be used to construct the similarity matrix. Based on the association testing of 6228 almost uncorrelated CVs from the GAW18 data, we find that PC-based models were capable of taking into account the sample correlations and worked as well as the LMM. This result is different from the claim made in Price et al  that PC-based models do not model family structure or cryptic relatedness. When using CVs to construct the similarity matrix, the top few PCs from the IBS matrix and the covariance matrix yield similar results. But when using RVs, the top few PCs from the IBS matrix are slightly better than those from the covariance matrix. LMM implemented by EMMAX is generally as effective as anticipated, although sometimes it can be conservative.
One limitation in our study was that in GAW18 data, there is no serious inflation in type I error for testing HTN1, even without any adjustment. Although our studies show a positive answer, more studies might be needed to confirm the effectiveness of a PC-based model for testing the binary trait.
We thank the GAW18 workshop for offering us such a complete data set to practice our methods. YZ and WP were supported by National Institutes of Health grants R21DK089351, R01HL65462, R01HL105397, and R01GM081535. 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.
- Patterson N, Price AL, Reich D: Population structure and eigenanalysis. PLoS Genet. 2006, 2: e190-10.1371/journal.pgen.0020190.PubMed CentralView ArticlePubMedGoogle Scholar
- Price AL, Zaitlen NA, Reich D, Patterson N: New approaches to population stratification in genome-wide association studies. Nat Rev Genet. 2010, 11: 459-463.PubMed CentralView ArticlePubMedGoogle Scholar
- Yu J, Pressoir G, Briggs WH, Bi IV, Yamasaki M, Doebley JF, McMullen MD, Gaut BS, Nielsen DM, Holland JB: A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat Genet. 2005, 38: 203-208.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
- Kang HM, Sui JH, Zaitlen NA, Kong S, 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
- Siu H, Lin L, Xiong M: Manifold learning for human population structure studies. PLoS One. 2012, 7: E29901-10.1371/journal.pone.0029901.PubMed CentralView ArticlePubMedGoogle Scholar
- Baye TM, He H, Ding L, Kurowski BG, Zhang X, Martin LJ: Population structure analysis using rare and common functional variants. BMC Proc. 2011, 5 (Suppl 9): S8-10.1186/1753-6561-5-S9-S8.PubMed CentralView ArticlePubMedGoogle Scholar
- Zhang Y, Guan W, Pan W: Adjustment for population stratification via principal components in association analysis of rare variants. Genet Epidemiol. 2013, 37: 99-109. 10.1002/gepi.21691.PubMed CentralView ArticlePubMedGoogle Scholar
- Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, Maller J, Sklar P, De Bakker PIW, Daly MJ: 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.