Statistical Power to Detect Genetic (Co)Variance of Complex Traits Using SNP Data in Unrelated Samples

1Genome-wide association studies (GWAS) have been extremely successfully in identifying genetic variants associated with complex traits and diseases in humans. In GWAS, hundreds of thousands or millions of SNPs are tested one by one for statistical evidence of association with a trait, and to avoid false positive discoveries due to the very large number of statistical tests being conducted, usually a very stringent p-value threshold, e.g. 5×10−8, is used to report a significant finding. Therefore, if there are many genes each with a small effect affecting the trait, most of these genetic variants will fail to pass the stringent threshold and remain undetected. This is one of the explanations of the ‘missing heritability’ question, that genetic variants identified from GWAS so far explain a fraction of the heritability for complex traits. We proposed a method, which is able to estimate the total amount of variance explained by all SNPs together without testing the SNPs individually for a quantitative trait, and subsequently extended it to the estimation of missing heritability for binary disease data from ascertained case-control studies. The analyses until recently only included common SNPs (e.g. minor allele frequency >0.01). The estimate quantifies the overall contribution from the additive effects of all SNPs, which is the upper limit of the proportion of variance that is captured by the additive effects of the set of SNPs used in the estimation, and is also the lower limit of the narrow-sense heritability of the trait. We also extended the method to estimate the genetic correlation between two traits using SNP data,. In contrast to the traditional (co)variance estimation methods that rely on pedigree information (family/twin studies), our method uses unrelated samples from a general population and the genetic (co)variance is estimated using a genetic relationship matrix (GRM) estimated from SNPs. The estimate of genetic variance using SNP data in unrelated individuals is free of confounding from common environment effects shared between close relatives that are difficult to model in family-based analyses, and is directly comparable to results from GWAS, because both are based on the same experimental design. For multiple trait analysis, the SNP-based approach allows the estimation of the genetic correlation between complex traits measured on different samples, _ENREF_8. This is important in particular for estimating the genetic correlation between diseases because multiple diseases are unlikely to co-segregate in sufficiently large pedigrees to allow estimation using traditional pedigree design. The SNP-based method has the flexibility of estimating the genetic correlation between any two diseases using completely independent case-control data. Other methods to estimate genetic parameters from individual-level or summary GWAS data have also been reported –.

2We previously named the SNP-based method mentioned GREML, as a complement to GBLUP where variance components are assumed known, and have been implemented them in the software tool GCTA. One outstanding question is the statistical power of detecting genetic variation using the population-based estimation method, for example how many samples are required to achieve estimates that are sufficiently accurate to detect genetic (co)variance of complex traits. In this paper, we derive the sampling variance of the estimate of genetic (co)variance by analytical derivations and verify our derivations by simulations under a range of scenarios. We also provide an online tool for power calculation.

3The methods of using SNP data to estimate genetic variance in unrelated individuals have been detailed elsewhere,. In brief, given GWAS data, we can model the phenotype aswhere y is an N×1 vector of phenotypes with N being the sample size, g is an N×1 vector with each of its elements being the total genetic effect of an individual captured by all SNPs, and e is an N×1 vector of residuals. We have and, where is the genetic variance captured by all SNPs, A is the genetic relationship matrix (GRM) estimated from SNPs, is the residuals variance and I is an identity matrix. The genetic relationships, also known as ‘genomic relationships’ or ‘genetic similarity relationships’, are referenced to the current population, and so can be negative as they are distributed about a mean of zero. Equation (1) is a typical mixed linear model with, in which the variance components can be estimated using a restricted maximum likelihood (REML) approach,. The proportion of variance explained by all SNPs (SNP heritability) is defined as.

4For power calculation, we need to know the sampling variance of the estimate of, i.e.. In practice, the asymptotic sampling variance (standard error squared) of a variance component is calculated from a diagonal element of the inverse of the information matrix in maximum likelihood analysis –. Each element of the information matrix, however, comprises complex forms of matrix algebra including a matrix inverse. It is therefore unfeasible to derive directly from the inverse of the information matrix. We show below an equivalent approach to obtain under the simple regression framework.

5For unrelated individuals, where the phenotypic correlation between individuals is small, mixed linear model analysis using the REML approach is asymptotically equivalent to simple regression analysis of pairwise phenotypic similarity/difference on pairwise genetic similarity, as measured by identity-by-descent (IBD) or identity-by-state (IBS) at genome-wide markers –. Under such circumstance, a regression of the cross-product of the phenotypes is equivalent to using both the squared difference and squared sum of the pairwise phenotypes, and using the cross-product is equivalent to using maximum likelihood. The model for the regression-based analysis can be written aswhere with and being the phenotypes of individuals i and j, is the ij-th element of the GRM A, and is the residual of this regression. There are observations (contrasts) in the regression. The regression coefficient b is equivalent to becauseIn such a simple regression, the sampling variance of the estimate of the regression coefficient isIf the samples are unrelated and the phenotypes have been standardized with mean of 0 and variance of 1, then and. Since is small, there is hardly any variance in that can be explained by so that. We therefore haveUnder circumstances when is large, for example when the GRM is calculated from pedigree data, a substantial proportion of variance in could be explained by, so that will be smaller than and the sampling variance of estimate of genetic variance will be reduced accordingly. In general, and the residual variance in equation (2) depend on the number of SNP that are used to calculate the GRM and their correlation structure. Although can be calculated empirically from the data, theoretical work suggest it is approximately 2×10−5 for genome-wide coverage of common SNPs in human populations. Since the phenotypic variance is usually estimated with very high precision,This suggests that the standard error (SE) of depends only on sample size, and is approximately. We show by simulations based on real genotype data (Text S1) that this approximation is very accurate (Figure 1 and Table S1). The SEs calculated from the approximation theory are also highly consistent with those reported from our previous studies for human height and body mass index (BMI). For example, the reported SE of for height was 0.083 using 3925 unrelated samples and 0.029 for both height and BMI, irrespective to, using 11586 unrelated samples, and the SE calculated from the approximation theory is 0.081 for N = 3925 and 0.027 for N = 11586.

6For a bivariate analysis where the two traits are measured on the same individuals, the mixed linear model can be written as where y 1 and y 2 are N×1 vector of phenotypes, g 1 and g 2 are N×1 vectors of genetic effects with and, e 1 and e 2 are N×1 vectors of residuals with and, and N is the sample size. The variance covariance matrix iswhere is the genetic covariance between the two traits and is the residual covariance. The genetic variance and covariance components can also be estimated using REML. The genetic correlation is estimated as. Since is a non-linear function of, and, there is no explicit derivation for. Reeve (1955) and Robertson (1959) provided an approximation of in the context of balanced pedigree design as, and Koots and Gibson (1996) proposed a modified version as, where is the phenotypic correlation between the traits. However, both approximations have an unsatisfying property that will approach 0 if or is close to 1. We derived an approximation, which does not have this problem (Text S2), i.e.As described above, for a GRM estimated from common SNPs in unrelated individuals in human populations, therefore. When, i.e. the traits are completely independent,. We tested equation (7) by simulations based on real genotype data (Text S1). The simulation results suggest that the approximation is reasonably accurate (Table 1). For real data analysis, we previously estimated the genetic correlation between intelligence at age 11 years and in old age of 0.62 with a SE of 0.23 using 1729 samples, consistent with the predicted SE of 0.22 from the approximation theory.

7For a bivariate analysis where the two traits are measured on different sets of individuals, e.g. height in males and blood pressure in females, the variance-covariance matrix is where y 1 is an N 1×1 vector of phenotypes in sample set #1 (e.g. males), and y 2 is an N 2×1 vector of phenotypes for in sample set #2 (e.g. females), with N 1 and N 2 being the sample sizes of the two sets. A 1 is an N 1×N 1 GRM for individuals in sample set #1, A 2 is an N 2×N 2 GRM in sample set #2 and A 12 is an N 1×N 2 GRM between the two sets of samples. and are the genetic variance for the two traits. and are the residual variances with the corresponding identify matrix I 1 and I 2. is the genetic covariance between traits. Since the two traits are measured on different sets of samples, the residual covariance is ignored because it is assumed that there is no covariance between the unrelated individuals apart from that caused by genetic factors. The genetic correlation is also estimated as, however, the sampling variance of is different from that described above. Since the traits are measured in different sets of samples,. Therefore, from a second order Taylor series approximation This approximation involves the sampling variance of. We show below an equivalent approach to obtain.

8Analogous to the univariate analysis, estimation of genetic covariance by a bivariate mixed linear model analysis is asymptotically equivalent to the following linear regression modelwhere i.e. the product of phenotypes between the i-th individual in set #1 and the j-th individual in set #2, and is the ij-th element of the GRM A 12, i.e. the genetic relationship between the i-th individual in set 1 and the j-th individual in sample set #2. The regression coefficient is equivalent to genetic covariance between the two traits becauseIf the two sample sets are independent and phenotypes for both traits have been standardized with mean of 0 and variance of 1, then and. Since is small,. We then have.

9We know from the derivations above that and. For unrelated individuals sampled from the same population,, we therefore getThis was also tested by simulations (Text S1) and the approximated standard errors were highly consistent with those observed from simulations, especially when sample size was large (Table 1). When, i.e. two traits are completely independent, for traits measured on the same sample, and for traits measured on different samples. Therefore, for independent traits, the ratio of sampling variance of genetic correlation between the two traits measured on the same sample to that on different samples is simply.

10For case-control studies, the proportion of variance in case-control status (0 or 1) that is explained by all SNPs on the observed scale can be estimated using a linear model. Therefore, the same approximations to the sampling variance of genetic variance and genetic correlation for quantitative traits can be applied directly to case-control studies. As shown in equation (5), in a univariate analysis, the sampling variance of SNP-based heritability depends only on sample size and variance in genetic relatedness, independent of the properties of the phenotype, so that var is also approximately in a case-control study with N being the total number of cases and controls. We show in Table 2 that the observed standard errors of the estimates of from published studies are highly consistent with those predicted from our approximation theory.

11To calculate power, however, we would need to specify, which is a parameter with non-intuitive properties, and depends on the prevalence of the disease in the population (K), the proportion of variance in disease liability that is captured by the SNPs at population level, and the proportional of cases in the sample (v). For this reason we define as the variance explained by all SNPs at the population level on the unobserved underlying scale of disease liability, and use a linear transformation to transform to on the liability scale, i.e. and with. We then getwhere N case is the number of cases, N control is the number of controls, and i is the selection intensity which is a function of K. We illustrate in Figure 2 the dependency of the SE of on disease prevalence (K) and proportion of cases in the sample (v) due to the transformation.

12As shown in equation (10), in a bivariate analysis where traits are measured on different sets of samples, the sampling variance of genetic correlation depends on sample sizes, trait heritabilities and the genetic correlation parameter, which is also independent of the properties of the phenotypes. Therefore, in a bivariate analysis of two independent case-control disease studies,where N 1 and N 2 are the total numbers of cases and controls of the two case-control studies, respectively. This also applies to a bivariate analysis of a quantitative trait and a cases-control disease study on different sets of samples, i.e.These two equations can also be expressed with respect to, given (see above). We show in Table 3 that the reported SEs of from bivariate analyses of psychiatric diseases are also highly in line with the predicted SEs from the approximation theory.

13Statistical power is calculated from the population value of the parameter and its sampling variance, which was derived above. If the parameter is θ, where θ is either the proportion of phenotypic variance captured by SNPs in the univariate case or the genetic correlation in the bivariate case, then is asymptotically distributed as a non-central χ 2 with 1 degree of freedom and non-centrality parameter (NCP) of. Given λ and the type-I error rate of α, statistical power is the probability that a non-central χ 2 variable is larger than the central χ 2 threshold that is determined by α. We show in Figure 3 the statistical power based on the sampling variance from our approximation theories to detect in a univariate case and in a bivariate case under a range of scenarios. For example, for a quantitative trait, approximately 8900, 4500, 3000 and 2300 independent individuals are required to detect of 0.1, 0.2, 0.3 and 0.4 with >80% power at a type-I error rate of 0.05, respectively. For two quantitative traits measured on the same sample, approximately 7000, 4700, 2500 and 1600 independent individuals are required to detect of 0.2, 0.4, 0.6 and 0.8 with >80% at a type-I error rate of 0.05, respectively.

14We have also developed an online calculator (GCTA Power Calculator, http://spark.rstudio.com/ctgg/gctaPower), as part of the GCTA software package (http://ctgg.qbi.uq.edu.au/software/gcta), using R-Shiny (http://shiny.rstudio.org) to calculate the SE of genetic variance or genetic correlation and statistical power given user-defined parameters.

15We have derived the approximate sampling variance of the estimate of variance explained by all common SNPs for a quantitative trait or case-control study of a disease, and genetic correlation between two quantitative traits, between two diseases, or between a quantitative trait and a disease, using genome-wide SNP data in unrelated individuals. We believe that the derivations and the online tool will be helpful for researchers to determine how many samples are required to detect (or ) and to estimate (or ) with adequate precision before collecting the genotype data.

16The sampling variance of for a complex trait is inversely proportional to sample size (N) and the variance in SNP-based genetic relatedness, and independent of. The sampling variance of between two complex traits is a function of, N of the two samples, of the two traits and when the traits are measured on different samples, and further depends on the phenotypic correlation when traits are measured on the same samples. All the approximation theories apply to case-control studies of diseases since the case-control data can be analysed using a linear model on the observed 0–1 scale. The sampling variance for the estimate on the observed scale can then be transformed to that on the underlying liability scale using well-established theory. The standard errors (square root of sampling variance) of either or observed in published studies were all highly consistently with those predicted from our approximation theories, which were also confirmed by simulations based on real genotype data.

17Analytical expressions for the sampling variance of the estimates of genetic (co)variance from pedigree analyses have been around for over 50 years,, and statistical power can be derived from these by using the sampling variance and population value of the parameter. However, these expressions are typically for specific structured pedigrees, such as fullsib or halfsib families or twin pairs. There are to our knowledge no simple approximations for general pedigrees, because the inverse of the variance-covariance matrix is required and this is conditional on the actual pedigree structure. The sampling variance of the estimated parameters in a general complex pedigree is usually derived post hoc after the analysis has been performed.

18Methods for calculating the power of detecting quantitative trait loci (QTL) in family-based linkage studies have been investigated extensively in the past two decades –,. These methods were developed to calculate the power of detecting a QTL but can be generalized for variance components estimation, e.g. estimating the genetic variance using pedigree information. The non-centrality parameter of the test-statistic from a maximum likelihood analysis of variance components is, where L is the likelihood function, and and are the variance covariance matrix under the null and alternative hypotheses respectively,. For a specific balanced pedigree design, e.g. fullsibs or nuclear families, the determinant (or inverse) of the V matrix can be computed explicitly, so that the NCP can be calculated without making approximation,. For an arbitrary pedigree, can be calculated approximately using Taylor expansions given the variance in family relatedness,. Therefore, all these methods explicitly or implicitly require a known pedigree. When the correlations between relatives are small, the first order approximation of the NCP in Rijsdijk et al can be written in our notations as, which is the same as we derived (i.e., see Equation (4) for ), even though our deviations were based on least squares regression analysis in unrelated samples whereas the derivations in Rijsdijk et al were based on maximum likelihood approach in family data. This approximation is reasonably accurate when correlations between relatives are small for a pedigree-based design, which is not an issue for a population-based design where the genetic relationships between unrelated samples are very small as demonstrated in Yang et al. We show by simulations (Text S1) that for a univariate analysis the LRT statistics calculated based on REML are highly consistent with the chi-squared test-statistics calculated by the Wald test using the sampling variance either observed from the simulations or predicted from our approximation theory (Figure S1).

19For a given population, a set of common SNPs and the method of calculating the genetic relationship matrix that we have used here, is a fixed quantity because it depends only on effective population size of the human populations. We used, which was calculated from theory based upon an effective population size of 10,000. Variance in genetic relatedness (and therefore power of detection) can decrease by including many rare SNPs in calculating the GRM because adding more rare SNPs increases the effective population size reflecting recent population expansion. The variance in relatedness can also increase by sampling closer relatives (see below for more discussion) or, for example, by creating a relationship matrix based upon haplotype information. Modifying the GRM can also affect the variance of the off-diagonal elements. For example by applying a weighting of SNPs depending on linkage disequilibrium the variance in the estimates of genetic relationships will decrease so that the sampling variance of the estimate of SNP-based heritability will be increased. Although we derive the theory and show the results based on the SNPs on the whole genome, our approximation theories are also applicable in analyses using a subset of SNPs, e.g. SNPs from a single chromosome. In that case, used in the approximation equations should be either observed empirically from data or derived from theory based on the subset of SNPs.

20If there are unknown related samples in the data (cryptic relatedness), will possibly be inflated due to shared environment between close relatives and/or the effects of causal variants in LD with the SNPs but captured by family relatedness, and will be deflated due to the increase of. In fact, the interpretation of changes if there is a substantial proportion of close relatives in the data,. This, however, affects GWAS result in a similar way, where the SE of the estimate of a SNP effect from a single SNP analysis (e.g. linear regression for a quantitative trait and logistic regression for a case-control study) will be deflated, causing an inflation of the test-statistics GWAS (often called “genomic inflation” ). For the estimation of using all common SNPs, to avoid possible confounding from shared environments and uncaptured causal variants, we suggested in Yang et al. (2010) a stringent threshold, i.e. 0.025, to remove cryptic relatedness from the data so that the estimate of can be compared directly to the results from GWAS in response to the “missing heritability” problem. In practice, observing a much smaller SE of using all common SNPs than that predicted from theory is a caveat suggesting substantial cryptic relatedness remaining in the data.

21Using the same experimental design of a sample of conventionally unrelated individuals, the experimenter can increase power by increasing sample size. Fortunately, power increases quadratically with sample size because every new sample is contrasted with all existing samples. The sampling variance of the estimate of the genetic correlation is generally much larger than that of the proportion of variance explained from a univariate analysis, consistent with the theory of the sampling variance of genetic correlations in pedigree designs.