A new method for estimating effect size distribution and heritability from genome-wide association summary results

Abstract

Accurately estimating the distribution and heritability of SNP effects across the genome could help explain the mystery of missing heritability. In this study, we propose a novel statistical method for estimating the distribution and heritability of SNP effects from genome-wide association studies (GWASs), and compare its performance to several existing methods using both simulations and real data. Specifically, we study the full range of GWAS summary results and link observed p values and unobserved effect sizes by (non-central) Chi-square distribution. By modeling the observed full set of association signals using a multinomial distribution, we build a likelihood function of SNP effect sizes using parametric and non-parametric maximum likelihood frameworks. Simulation studies show that the proposed method can accurately estimate effect sizes and the number of associated SNPs. As real applications, we analyze publicly available GWAS summary results for height, body mass index (BMI), and bone mineral density (BMD). Our analyses show that there are over 10,000 SNPs that might be associated with height, and the total heritability attributable to these SNPs exceeds 70 %. The heritabilities for BMI and BMD are ~10 and ~15 %, respectively. The results indicate that the proposed method has the potential to improve the accuracy of estimates of heritability and effect size for common SNPs in large-scale GWAS meta-analyses. These improved estimates may contribute to an enhanced understanding of the genetic basis of complex traits.

Appendices

Appendix 1: The form of Chi-squared non-central parameter (NCP) as a function of effect size (ES)

Under large sample size, the Wald test of association is approximately equivalent to the likelihood ratio test (LRT). Below, we will derive the form of Chi-squared NCP value λ with the LRT instead.

For a particular SNP site, denote the two alleles as A and a, and their frequencies as f and 1 − f. The three genotypes, AA, Aa, and aa are encoded as 0, 1, and 2 under the additive mode of inheritance. Under the Hardy–Weinberg equilibrium, their frequencies are f ², 2f(1 − f) and (1 − f)².

The linear regression model to test the association between the phenotype Y and the genotype G is

$$y = \alpha + \beta \cdot g + \varepsilon,$$

where α is the regression intercept, β is regression coefficient, ε is the normally distributed random error ε ~ N(0, σ ²). The parameters of the above model are θ = (α, β, σ ²), and the hypothesis to be tested is β = 0.

Under the alternative hypothesis, suppose the true model parameters are θ* = (α*, β*, σ*²). The distribution of the phenotype Y is dependent upon G in that y ~ N(α* + β*g, σ*²). The effect size e of the SNP is measured by the proportion of variation in Y explained by G,

$$e = \frac{{\text{var} (\beta^{*} g)}}{{\text{var} (\beta^{*} g) + \text{var} (e)}} = \frac{{\beta^{*2} \cdot 2f(1 - f)}}{{\beta^{*2} \cdot 2f(1 - f) + \sigma^{*2}}}$$

To calculate the NCP, we refer to the reference (Gauderman 2003) and take Y as random variable. Specifically, we calculate expected log-likelihood functions with respect to Y and G under the alternative and null hypotheses, respectively. Subtraction of the two maximized expected log-likelihoods composed the element of NCP.

For a single subject, the expected log-likelihood function l is defined as

$$l = \sum _{G} \smallint _{Y} { \ln } \left\{{L\left({{\varvec{\uptheta}}; y, g} \right)} \right\}p\left({y, g|{\varvec{\uptheta}}^{*}} \right)dY,$$

where $L({\varvec{\uptheta}};y,g) \propto p(y|g,{\varvec{\uptheta}}) = \frac{1}{{\sqrt {2\pi \sigma^{2}}}}\exp \left\{{- \frac{{(y - \alpha - \beta g)^{2}}}{{2\sigma^{2}}}} \right\}$, and,

$$p(y,g|{\varvec{\uptheta}}^{*}) \propto p(y|g,{\varvec{\uptheta}}^{*})p(g|{\varvec{\uptheta}}^{*}) = \frac{1}{{\sqrt {2\pi \sigma^{*2}}}}\exp \left\{{- \frac{{(y - \alpha^{*} - \beta^{*} g)^{2}}}{{2\sigma^{*2}}}} \right\}p(g)$$

Under the alternative hypothesis,

$$\begin{aligned} l_{1} & = \sum\limits_{G} {\int\limits_{Y} {\ln \left[{\frac{1}{{\sqrt {2\pi \sigma^{2}}}}\exp \left\{{- \frac{{(y - \alpha - \beta g)^{2}}}{{2\sigma^{2}}}} \right\}} \right]\frac{1}{{\sqrt {2\pi \sigma^{*2}}}}\exp \left\{{- \frac{{(y - \alpha^{*} - \beta^{*} g)^{2}}}{{2\sigma^{*2}}}} \right\}p(g)dY}} \\ &= \sum\limits_{G} {\int\limits_{Y} {\left\{{- \frac{1}{2}\ln (2\pi \sigma^{2}) - \frac{{(y - \alpha - \beta g)^{2}}}{{2\sigma^{2}}}} \right\}\frac{1}{{\sqrt {2\pi \sigma^{*2}}}}\exp \left\{{- \frac{{(y - \alpha^{*} - \beta^{*} g)^{2}}}{{2\sigma^{*2}}}} \right\}p(g)dY}} \\ &= - \frac{1}{2}\ln (2\pi \sigma^{2}) - \frac{{f^{2} (\alpha^{*} - \alpha)^{2} + 2f(1 - f)(\alpha^{*} + \beta^{*} - \alpha - \beta)^{2} + (1 - f)^{2} (\alpha^{*} + 2\beta^{*} - \alpha - 2\beta)^{2} + \sigma^{*2}}}{{2\sigma^{2}}} \\ \end{aligned}$$

Maximizing l ₁ by partial derivatives with respect to α, β, and σ ² gets the MLE estimates,

$$\hat{\alpha}_{1} = \alpha^{*},\hat{\beta}_{1} = \beta^{*},\hat{\sigma}_{1}^{2} = \sigma^{*2},\quad {\text{and}}\quad \hat{l}_{1} = - \frac{1}{2}\ln (2\pi \sigma^{*2}) - \frac{1}{2}$$

Under the null hypothesis, β is restricted to be zero,

$$\begin{aligned} l_{0} & = \sum\limits_{G} {\int\limits_{Y} {\ln \left[{\frac{1}{{\sqrt {2\pi \sigma^{2}}}}\exp \left\{{- \frac{{(y - \alpha)^{2}}}{{2\sigma^{2}}}} \right\}} \right]\frac{1}{{\sqrt {2\pi \sigma^{*2}}}}\exp \left\{{- \frac{{(y - \alpha^{*} - \beta^{*} g)^{2}}}{{2\sigma^{*2}}}} \right\}p(g)dY}} \\ & = \sum\limits_{G} {\int\limits_{Y} {\left\{{- \frac{1}{2}\ln (2\pi \sigma^{2}) - \frac{{(y - \alpha)^{2}}}{{2\sigma^{2}}}} \right\}\frac{1}{{\sqrt {2\pi \sigma^{*2}}}}\exp \left\{{- \frac{{(y - \alpha^{*} - \beta^{*} g)^{2}}}{{2\sigma^{*2}}}} \right\}p(g)dY}} \\ & = - \frac{1}{2}\ln (2\pi \sigma^{2}) - \frac{{f^{2} (\alpha^{*} - \alpha)^{2} + 2f(1 - f)(\alpha^{*} + \beta^{*} - \alpha)^{2} + (1 - f)^{2} (\alpha^{*} + 2\beta^{*} - \alpha)^{2} + \sigma^{*2}}}{{2\sigma^{2}}} \\ \end{aligned}$$

Maximizing l ₀ by partial derivatives with respect to α and σ ² gets the MLE estimates,

$$\hat{\alpha}_{0} = \alpha^{*} + 2\beta^{*} (1 - f),\hat{\sigma}_{0}^{2} = \sigma^{*2} + 2f(1 - f)\beta^{*2}$$

Recalling the expression of e gets

$$\hat{\sigma}_{0}^{2} = \sigma^{*2} \left({\frac{1}{1 - e}} \right),\quad {\text{and}}$$

$$\hat{l}_{0} = - \frac{1}{2}\ln (2\pi \sigma^{*2}) - \frac{1}{2}\ln \left({\frac{1}{1 - e}} \right) - \frac{1}{2}$$

Therefore, for a sample of N subjects,

$$\lambda = 2 N(\hat{l}_{1} - \hat{l}_{0}) = N\ln \left({\frac{1}{1 - e}} \right).$$

Appendix 2: Integration form for gamma distribution

For gamma distribution family gamma(shape = k, scale = θ), the integration term has the following form

$$\begin{aligned} \int_{0}^{1} {f_{\lambda} (q_{r})p_{{\varvec{\uptheta}}} (e)de} & = \sum\limits_{i = 0}^{\infty} {\frac{{q^{{\frac{1 + 2i}{2} - 1}} \cdot e^{{- \frac{q}{2}}} \cdot \left({\frac{N}{2}} \right)^{i} \cdot \left({\frac{1}{{\frac{N}{2} + \frac{1}{\theta}}}} \right)^{i + k} \cdot \varGamma (i + k)}}{{i! \cdot 2^{{\frac{1 + 2i}{2}}} \cdot \varGamma \left({\frac{1 + 2i}{2}} \right) \cdot \theta^{k} \cdot \varGamma (k)}}} \\ & = \sum\limits_{i = 0}^{\infty} {\exp \left\{{\left({\frac{1 + 2i}{2} - 1} \right)\log (q) - \frac{q}{2} + i\log \left({\frac{N}{2}} \right) + (i + k)\log \left({\frac{2\theta}{2 + N\theta}} \right) + \log (\varGamma (i + k)) - \log (i!) - \left({\frac{1 + 2i}{2}} \right)\log (2) - \log \left({\varGamma \left({\frac{1 + 2i}{2}} \right)} \right) - k\log (\theta) - \log (\varGamma (k))} \right\}} \\ \end{aligned}.$$

Each term in the above equation is easily to compute. The summation will start from i = 0 and will usually converge after a small number of iterations (small i).