Beyond SNP heritability: Polygenicity and discoverability of phenotypes estimated with a univariate Gaussian mixture model

Summary

Estimating the polygenicity (proportion of causally associated single nucleotide polymorphisms (SNPs)) and discoverability (effect size variance) of causal SNPs for human traits is currently of considerable interest. SNP-heritability is proportional to the product of these quantities. We present a basic model, using detailed linkage disequilibrium structure from a reference panel of 11 million SNPs, to estimate these quantities from genome-wide association studies (GWAS) summary statistics. We apply the model to diverse phenotypes and validate the implementation with simulations. We find model polygenicities (as a fraction of the reference panel) ranging from ≃ 2 × 10⁻⁵ to ≃ 4 × 10⁻³, with discoverabilities similarly ranging over two orders of magnitude. A power analysis allows us to estimate the proportions of phenotypic variance explained additively by causal SNPs reaching genome-wide significance at current sample sizes, and map out sample sizes required to explain larger portions of additive SNP heritability. The model also allows for estimating residual inflation (or deflation from over-correcting of z-scores), and assessing compatibility of replication and discovery GWAS summary statistics.

Keywords:

Heredity – Molecular genetics – Heterozygosity – Genome-wide association studies – Alzheimer's disease – Schizophrenia – Amyotrophic lateral sclerosis – Fourier analysis

Introduction

The genetic components of complex human traits and diseases arise from hundreds to likely many thousands of single nucleotide polymorphisms (SNPs) [1], most of which have weak effects. As sample sizes increase, more of the associated SNPs are identifiable (they reach genome-wide significance), though power for discovery varies widely across phenotypes. Of particular interest are estimating the proportion of common SNPs from a reference panel (polygenicity) involved in any particular phenotype; their effective strength of association (discoverability, or causal effect size variance); the proportion of variation in susceptibility, or phenotypic variation, captured additively by all common causal SNPs (approximately, the narrow sense heritability), and the fraction of that captured by genome-wide significant SNPs—all of which are active areas of research [2–9]. The effects of population structure [10], combined with high polygenicity and linkage disequilibrium (LD), leading to spurious degrees of SNP association, or inflation, considerably complicate matters, and are also areas of much focus [11–13]. Despite these challenges, there have been recent significant advances in the development of mathematical models of polygenic architecture based on GWAS [14, 15]. One of the advantages of these models is that they can be used for power estimation in human phenotypes, enabling prediction of the capabilities of future GWAS.

Here, in a unified approach explicitly taking into account LD, we present a model relying on genome-wide association studies (GWAS) summary statistics (z-scores for SNP associations with a phenotype [16]) to estimate polygenicity (π₁, the proportion of causal variants in the underlying reference panel of approximately 11 million SNPs from a sample size of 503) and discoverability (σ β 2, the causal effect size variance), as well as elevation of z-scores due to any residual inflation of the z-scores arising from variance distortion (σ 0 2, which for example can be induced by cryptic relatedness), which remains a concern in large-scale studies [10]. We estimate π₁, σ β 2, and σ 0 2, by postulating a z-score probability distribution function (pdf) that explicitly depends on them, and fitting it to the actual distribution of GWAS z-scores.

Estimates of polygenicity and discoverability allow one to estimate compound quantities, like narrow-sense heritability captured by the SNPs [17]; to predict the power of larger-scale GWAS to discover genome-wide significant loci; and to understand why some phenotypes have higher power for SNP discovery and proportion of heritability explained than other phenotypes.

In previous work [18] we presented a related model that treated the overall effects of LD on z-scores in an approximate way. Here we take the details of LD explicitly into consideration, resulting in a conceptually more basic model to predict the distribution of z-scores. We apply the model to multiple phenotype datasets, in each case estimating the three model parameters and auxiliary quantities, including the overall inflation factor λ, (traditionally referred to as genomic control [19]) and narrow sense heritability, h². We also perform extensive simulations on genotypes with realistic LD structure in order to validate the interpretation of the model parameters. A discussion of the relation of the present paper to other work is provided in the first section of the S1 Appendix (pp. S2-S3).

Materials and methods

Overview

Our basic model is a simple postulate for the distribution of causal effects (denoted β below) [20]. Our model assumes that only a fraction of all SNPs are in some sense causally related to any given phenotype. We work with a reference panel of approximately 11 million SNPs with 503 samples, and assume that all common causal SNPs (minor allele frequency (MAF) > 0.002) are contained in it. Any given GWAS will have z-scores for a subset of these reference SNPs (we use the term “typed” below to refer to GWAS SNPs with z-scores, whether they were directly genotyped or their genotype was imputed). When a z-score partially involves a latent causal component (i.e., not pure noise), we assume that it arises through LD with neighboring causal SNPs, or that it itself is causal.

We construct a pdf for z-scores that directly follows from the underlying distribution of effects. For any given typed SNP’s z-score, it is dependent on the other SNPs the focal SNP is in LD with (SNPs that are “tagged” by the focal SNP), taking into account their LD with the focal SNP and their heterozygosity (i.e., it depends not just on the focal typed SNP’s total LD and heterozygosity, but also on the distribution of neighboring reference SNPs in LD with it and their heterozygosities). We present two ways of constructing the model pdf for z-scores, using multinomial expansion, and using convolution. The former is perhaps more intuitive, but the latter is more numerically tractable, yielding an exact solution, and is used here to obtain all reported results. The problem then is finding the three model parameters that give a maximum likelihood best fit for the model’s prediction of the distribution of z-scores to the actual distribution of z-scores. Because we are fitting three parameters typically using ≳10⁶ data points, it is appropriate to incorporate some data reduction to facilitate the computations. To that end, we bin the data (z-scores) into a 10 × 10 grid of heterozygosity-by-total LD (having tested different grid sizes to ensure convergence of results). Also, when building the LD and heterozygosity structures of reference SNPs, we fine-grained the LD range (0 ≤ r² ≤ 1), again ensuring that bins were small enough that results were well converged. To fit the model to the data we bin the z-scores (within each heterozygosity/total LD window) and calculate the multinomial probability for having the actual distribution of z-scores (numbers of z-scores in the z-score bins) given the model pdf for the distribution of z-scores, and adjusting the model parameters using a multidimensional unconstrained nonlinear minimization (Nelder-Mead), so as to maximize the likelihood of the data, given the parameters.

A visual summary of the predicted and actual distribution of z-scores is obtained by making quantile-quantile plots showing, for a wide range of significance thresholds going well beyond genome-wide significance, the proportion (x-axis) of typed SNPs exceeding any given threshold (y-axis) in the range. It is important also to assess the quantile-quantile sub-plots for SNPs in the heterozygosity-by-total LD grid elements (see S1 Appendix).

With the pdf in hand, various quantities can be calculated: the number of causal SNPs; the expected genetic effect (denoted δ below, where δ² is the non-centrality parameter of a Chi-squared distribution) at the current sample size for a typed SNP given the SNP’s z-score and its full LD and heterozygosity structure; the estimated SNP heritability, h S N P 2 (excluding contributions from rare reference SNPs, i.e., with MAF<0.2%); and the sample size required to explain any percentage of that with genome-wide significant SNPs. The model can easily be extended using a more complex distribution for the underlying β’s, with multiple-component mixtures for small and large effects, and incorporating selection pressure through both heterozygosity dependence on effect sizes and linkage disequilibrium dependence on the prior probability of a SNP’s being causal—issues we will address in future work.

The model: Probability distribution for z-scores

To establish notation, we consider a bi-allelic genetic variant, i, and let β_i denote the effect size of allele substitution of that variant on a given quantitative trait. We assume a simple additive generative model (simple linear regression, ignoring covariates) relating genotype to phenotype [18, 21]. That is, assume a linear vector equation (no summation over repeated indices)

for phenotype vector y over N samples (mean-centered and normalized to unit variance), mean-centered genotype vector g_i for the i^th of n SNPs (vector over samples of the additively coded number of reference alleles for the i^th variant), true fixed effect β_i (regression coefficient) for the SNP, and residual vector e_i containing the effects of all the other causal SNPs, the independent random environmental component, and random error. Variants with non-zero fixed effect β_i are said to be “causal”. For SNP i, the estimated simple linear regression coefficient is

where T denotes transpose and g i T g i / N = var ( g i ) = H i is the SNP’s heterozygosity (frequency of the heterozygous genotype): H_i = 2p_i(1 − p_i) where p_i is the frequency of either of the SNP’s alleles.

Consistent with the work of others [11, 15], we assume the causal SNPs are distributed randomly throughout the genome (an assumption that can be relaxed when explicitly considering different SNP categories, but that in the main is consistent with the additive variation explained by a given part of the genome being proportional to the length of DNA [22]). In a Bayesian approach, we assume that the parameter β for a SNP has a distribution (in that specific sense, this is similar to a random effects model), representing subjective information on β, not a distribution across tangible populations [23]. Specifically, we posit a normal distribution for β with variance given by a constant, σ β 2:

This is also how the β are distributed across the set of causal SNPs. Therefore, taking into account all SNPs (the remaining ones are all null by definition), this is equivalent to the two-component Gaussian mixture model we originally proposed [20]

where N ( 0 , 0 ) is the Dirac delta function, so that considering all SNPs, the net variance is var ( β ) = π 1 σ β 2. If there is no LD (and assuming no source of spurious inflation), the association z-score for a SNP with heterozygosity H can be decomposed into a fixed effect δ and a residual random environment and error term, ϵ∼ N ( 0 , 1 ), which is assumed to be independent of δ [18]:

with

so that

where

By construction, under null, i.e., when there is no genetic effect, δ = 0, so that var(ϵ) = 1.

If there is no source of variance distortion in the sample, but there is a source of bias in the summary statistics for a subset of markers (e.g., the sample is composed of two or more subpopulations with different allele frequencies for a subset of markers—pure population stratification in the sample [24]), the marginal distribution of an individual’s genotype at any of those markers will be inflated. The squared z-score for such a marker will then follow a non-central Chi-square distribution (with one degree of freedom); the non-centrality parameter will contain the causal genetic effect, if any, but biased up or down (confounding or loss of power, depending on the relative sign of the genetic effect and the SNP-specific bias term). The effect of bias shifts, arising for example due to stratification, is nontrivial, and currently not explicitly in our model; it is usually accounted for using standard methods [25].

Variance distortion in the distribution of z-scores can arise from cryptic relatedness in the sample (drawn from a population mixture with at least one subpopulation with identical-by-descent marker alleles, but no population stratification) [19]. If z_u denotes the uninflated z-scores, then the inflated z-scores are

where σ₀ ≥ 1 characterizes the inflation. Thus, from Eq 7, in the presence of inflation in the form of variance distortion

where σ ˜ β 2 ≡ σ 0 2 σ β 2, so that var ( δ ) = σ ˜ 2 ≡ σ ˜ β 2 N H and ∊ ∼ N ( 0 , σ 0 2 ). In the presence of variance distortion one is dealing with inflated random variables β ˜ ∼ N ( 0 , σ ˜ β 2 ), but we will drop the tilde on the β’s in what follows.

Since variance distortion leads to scaled z-scores [19], then, allowing for this effect in some of the extremely large data sets, we can assess the ability of the model to detect this inflation by artificially inflating the z-scores (Eq 9), and checking that the inflated σ ^ 0 2 is estimated correctly while the other parameter estimates remain unchanged.

Implicit in Eq 8 is approximating the denominator, 1 − q², of the χ² statistic non-centrality parameter to be 1, where q² is the proportion of phenotypic variance explained by the causal variant, i.e., q ≡ H β. So a more correct δ is

Taylor expanding in q and then taking the variance gives

The additional terms will be vanishingly small and so do not contribute in a distributional sense; (quasi-) Mendelian or outlier genetic effects represent an extreme scenario where the model is not expected to be accurate, but SNPs for such traits are by definition easily detectable. So Eq 8 remains valid for the polygenicity of complex traits.

Now consider the effects of LD on z-scores. The simple linear regression coefficient estimate for typed SNP i, β ^ i, and hence the GWAS z-score, implicitly incorporates contributions due to LD with neighboring causal SNPs. (A typed SNP is a SNP with a z-score, imputed or otherwise; generally these will compose a smaller set than that available in reference panels like 1000 Genomes used here for calculating the LD structure of typed SNPs.) In Eq 1, e_i = ∑_j≠i g_j β_j + ε, where g_j is the genotype vector for SNP j, β_j is its true regression coefficient, and ε is the independent true environmental and error residual vector (over the N samples). Thus, explicitly including all causal true β’s, Eq 2 becomes

(the sum over j now includes SNP i itself). This is the simple linear regression expansion of the estimated regression coefficient for SNP i in terms of the independent latent (true) causal effects and the latent environmental (plus error) component; β i ′ is the effective simple linear regression expression for the true genetic effect of SNP i, with contributions from neighboring causal SNPs mediated by LD. Note that g i T g j / N is simply cov(g_i, g_j), the covariance between genotypes for SNPs i and j. Since correlation is covariance normalized by the variances, β i ′ in Eq 13 can be written as

where r_ij is the correlation between genotypes at reference SNP j and typed SNP i. Then, from Eq 5, the z-score for the typed SNP’s association with the phenotype is given by:

We noted that in the absence of LD, the distribution of the residual in Eq 5 is assumed to be univariate normal. But in the presence of LD (Eq 15) there are induced correlations. Letting ϵ denote the vector of residuals (with element ϵ_i for SNP i, i = 1…n), and M denote the (sparse) n × n LD-r² matrix, then, ignoring inflation, ϵ ∼ N ( 0 , M ) [26]. Since the genotypes of two unrelated individuals are marginally independent, this multivariate normal distribution for ϵ is contingent on the summary statistics for all SNPs being determined from the same set of individuals, which generally is overwhelmingly, if not in fact entirely, the case (in the extreme, with an independent set of individuals for each SNP, M would be reduced to the identity matrix). A limitation of the present work is that we do not consider this complexity. This may account for the relatively minor misfit in the simulation results for cases of high polygenicity—see below.

Thus, for example, if the SNP itself is not causal but is in LD with k causal SNPs that all have heterozygosity H, and where its LD with each of these is the same, given by some value r² (0 < r² ≤ 1), then σ ˜ 2 in Eq 10 will be given by

For this idealized case, the marginal distribution, or pdf, of z-scores for a set of such associated SNPs is

where ϕ(⋅;μ, σ²) is the normal distribution with mean μ and variance σ², and H is shorthand for the LD and heterozygosity structure of such SNPs (in this case, denoting exactly k causal SNPs with LD given by r² and heterozygosity given by H). If a proportion α of all typed SNPs are similarly associated with the phenotype while the remaining proportion are all null (not causal and not in LD with causal SNPs), then the marginal distribution for all SNP z-scores is the Gaussian mixture

dropping the parameters for convenience.

For real genotypes, however, the LD and heterozygosity structure is far more complicated, and of course the causal SNPs are generally numerous and unknown. Thus, more generally, for each typed SNP H will be a two-dimensional histogram over LD (r²) and heterozygosity (H), each grid element giving the number of SNPs falling within the edges of that (r², H) bin. Alternatively, for each typed SNP it can be built as two one-dimensional histograms, one giving the LD structure (counts of neighboring SNPs in each LD r² bin), and the other giving, for each r² bin, the mean heterozygosity for those neighboring SNPs, which will be accurate for sufficiently fine binning—within a bin, the heterozygosities of the tagged referene SNPs wll be in a vary narrow range. We use the latter in what follows. We present two consistent ways of expressing the a posteriori pdf for z-scores, based on multinomial expansion and on convolution, that provide complementary views. The multinomial approach perhaps gives a more intuitive feel for the problem, but the convolution approach is considerably more tractable numerically and is used here to obtain all reporter results. All code used in the analyses, including simulations, is publicly available on GitHub [27].

Model PDF: Multinomial expansion

As in our previous work, we incorporate the model parameter π₁ for the fraction of all SNPs that are causal [18]. Additionally, we calculate the actual LD and heterozygosity structure for each SNP. That is, for each SNP we build a histogram of the numbers of other SNPs in LD with it for w equally-spaced r²-windows between r m i n 2 and 1 where r m i n 2 = 0 . 05 (approximately the noise floor for correlation when LD is calculated from the 503 samples in 1000 Genomes), and record the mean heterozygosity for each bin; as noted above, we use H as shorthand to represent all this. We find that w ≃ 20 is sufficient for converged results. For any given SNP, the set of SNPs thus determined to be in LD with it constitute its LD block, with their number given by n (LD with self is always 1, so n is at least 1). The pdf for z-scores, given N , H, and the three model parameters π₁, σ_β, σ₀, will then be given by the sum of Gaussians that are generalizations of Eq 17 for different combinations of numbers of causal SNPs among the w LD windows, each Gaussian scaled by the probability of the corresponding combination of causal SNPs among the LD windows, i.e., by the appropriate multinomial distribution term.

For w r²-windows, we must consider the possibilities where the typed SNP is in LD with all possible numbers of causal SNPs in each of these windows, or any combination thereof. There are thus w + 1 categories of SNPs: null SNPs (which r²-windows they are in is irrelevant), and causal SNPs, where it does matter which r²-windows they reside in. If window i has n_i SNPs (∑ i = 1 w n i = n) and mean heterozygosity H_i, and the overall fraction of SNPs that are causal is π₁, then the probability of having simultaneously k₀ null SNPs, k₁ causal SNPs in window 1, and so on through k_w causal SNPs in window w, for a nominal total of K causal SNPs (∑ i = 1 w k i = K and k₀ = n − K), is given by the multinomial distribution, which we denote M(k₀, …, k_w; n₀, …, n_w; π₁). For an LD block of n SNPs, the prior probability, p_i, for a SNP to be causal and in window i is the product of the independent prior probabilities of a SNP being causal and being in window i: p_i = π₁ n_i/n. The prior probability of being null (regardless of r²-window) is simply p₀ = (1 − π₁). The probability of a given breakdown k₀, …, k_w of the neighboring SNPs into the w + 1 categories is then given by

and the corresponding Gaussian is

For a SNP with LD and heterozygosity structure H, the pdf for its z-score, given N and the model parameters, is then given by summing over all possible numbers of total causal SNPs in LD with the SNP, and all possible distributions of those causal SNPs among the w r²-windows:

where K_max is bounded above by n. Note again that H is shorthand for the heterozygosity and linkage-disequilibrium structure of the SNP, giving the set {n_i} (as well as {H_i}), and hence, for a given π₁, p_i. Also there is the constraint ∑ i = 1 w k i = K on the second summation, and, for all i, max(k_i) = max(K, n_i), though generally K_max ≪ n_i. The number of ways of dividing K causal SNPs amongst w LD windows is given by the binomial coefficient ( a b ), where a ≡ K + w − 1 and b ≡ w − 1, so the number of terms in the second summation grows rapidly with K and w. However, because π₁ is small (often ≤10⁻³), the upper bound on the first summation over total number of potential causal SNPs K in the LD block for the SNP can be limited to K_max < min(20, n), even for large blocks with n ≃ 10³. That is,

Still, the number of terms is large; e.g., for K = 10 and w = 10 there are 92,378 terms.

For any given typed SNP (whose z-score we are trying to predict), it is important to emphasize that the specific LD r² and the heterozygosity of each underlying causal (reference) SNP tagged by it need to be taken into account, at least in an approximate sense that can be controlled to allow for arbitrary finessing giving converged results. This is the purpose of our w = 20 LD-r² windows, which inevitably leads to the multinomial expansion. Which window the causal SNP is in matters, leading to w+ 1 SNP categories, as noted above. Setting w = 1 would result in only a very rough approximation for the model pdf, reducing our multinomial to a binomial involving just two categories of SNPs: null and causal, with all causal SNPs treated the same, regardless of their LD with the tag SNP and their heterozygosity, as is done for the “M2” and “M3” models in [14]. The effects of this are demonstrated in the S1 Appendix (pp. S2-S3, and p. 14).

Model PDF: Convolution

From Eq 15, there exists an efficient procedure that allows for accurate calculation of a z-score’s a posteriori pdf (given the SNP’s heterozygosity and LD structure, and the phenotype’s model parameters). Any GWAS z-score is a sum of unobserved random variables (LD-mediated contributions from neighboring causal SNPs, and the additive environmental component), and the pdf for such a composite random variable is given by the convolution of the pdfs for the component random variables. Since convolution is associative, and the Fourier transform of the convolution of two functions is just the product of the individual Fourier transforms of the two functions, one can obtain the a posteriori pdf for z-scores as the inverse Fourier transform of the product of the Fourier transforms of the individual random variable components.

From Eq 15 z is a sum of correlation- and heterozygosity-weighted random variables {β_j} and the random variable ϵ, where {β_j} denotes the set of true causal parameters for each of the SNPs in LD with the typed SNP whose z-score is under consideration. The Fourier transform F(k) of a Gaussian f(x) = c × exp(−ax²) is F ( k ) = c π / a × exp ( - π 2 k 2 / a ). From Eq 4, for each SNP j in LD with the typed SNP (1 ≤ j ≤ b, where b is the typed SNP’s block size),

The Fourier transform (with variable k—see below) of the first term on the right hand side is

while that of the second term is simply (1 − π₁). Additionally, the environmental term is ϵ ∼ N ( 0 , σ 0 2 ) (ignoring LD-induced correlation, as noted earlier), and its Fourier transform is exp ( - 2 π 2 σ 0 2 k 2 ). For each typed SNP, one could construct the a posteriori pdf based on these Fourier transforms. However, it is more practical to use a coarse-grained representation of the data. Thus, in order to fit the model to a data set, we bin the typed SNPs whose z-scores comprise the data set into a two-dimensional heterozygosity/total LD grid (whose elements we denote “H-L” bins), and fit the model with respect to this coarse grid instead of with respect to every individual typed SNP z-score; in the section “Parameter Estimation” below we describe using a 10 × 10 grid. Additionally, for each H-L bin the LD r² and heterozygosity histogram structure for each typed SNP is built, using w_max equally-spaced r² bins for r m i n 2 ≤ r 2 ≤ 1 (this is a change in notation from the previous section: w_max here plays the role of w there; in what follows, w will be used as a running index, 1 ≤ w ≤ w_max); w_max = 20 is large enough to allow for converged results; r m i n 2 = 0 . 05 is generally small enough to capture true causal associations in weak LD while large enough to exclude spurious contributions to the pdf arising from estimates of r² that are non-zero due to noise. This points up a minor limitation of the model stemming from the small reference sample size (N_R = 503 for 1000 Genomes) from which H is built. Larger N_R would allow for more precision in handling very low LD (r² < 0.05), but this is an issue only for situations with extremely large σ β 2 (high heritability with low polygenicity) that we do not encounter for the 16 phenotypes we analyze here. In any case, this can be calibrated for using simulations.

We emphasize again that setting w_max = 1 would result in only an approximation for the model pdf (see “Relation to Other Work” in the S1 Appendix).

For any H-L bin with mean heterozygosity H and mean total LD L there will be an average LD and heterozygosity structure with a mean breakdown for the typed SNPs having n_w reference SNPs (not all of which necessarily are typed SNPs, i.e., have a z-score) with LD r² in the w th r² bin whose average heterozygosity is H_w. Thus, one can re-express z-scores for an H-L bin as

where β_j and ϵ are unobserved random variables.

In the spirit of the discrete Fourier transform (DFT), discretize the set of possible z-scores into the ordered set of n (equal to a power of 2) values z₁, …, z_n with equal spacing between neighbors given by Δz (z_n = −z₁ − Δz, and z_n/2+1 = 0). Taking z₁ = −38 allows for the minimum p-values of 5.8 × 10⁻³¹⁶ (near the numerical limit); with n = 2¹⁰, Δz = 0.0742. Given Δz, the Nyquist critical frequency is f c = 1 2 Δ z, so we consider the Fourier transform function for the z-score pdf at n discrete values k₁, …, k_n, with equal spacing between neighbors given by Δk, where k₁ = −f_c (k_n = −k₁ − Δk, and k_{n/2+ 1} = 0; the DFT pair Δz and Δk are related by ΔzΔk = 1/n). Define

(see Eq 24). Then the product (over r² bins) of Fourier transforms for the genetic contribution to z-scores, denoted G_j ≡ G(k_j), is

Recall that H denotes the LD and heterozygosity structure of a particular SNP (or representative SNP in an average sense for an H-L grid element), a shorthand for the set of values {n_w, H_w, L_w: w = 1, …, w_max} that characterize the SNP. Let M denote the set of model parameters. The Fourier transform of the environmental contribution, denoted E_j ≡ E(k_j), is

Let F_z = (G₁ E₁,…,G_nE_n) denote the vector of products of Fourier transform values, and let F - 1 denote the inverse Fourier transform operator. Then for the SNP in question, the vector of pdf values, pdf_z, for the uniformly discretized possible z-score outcomes z₁, …, z_n described above, i.e., pdf_z = (f₁,…,f_n) where f i ≡ pdf ( z i | H , M , N ), is

Thus, the i^th element pdfz i = f i is the a posteriori probability of obtaining a z-score value z_i for the SNP, given the SNP’s LD and heterozygosity structure, the model parameters, and the sample size.

Data preparation

For real phenotypes, we calculated SNP minor allele frequency (MAF) and LD between SNPs using the 1000 Genomes phase 3 data set for 503 subjects/samples of European ancestry [28–30]. In order to carry out realistic simulations (i.e., with realistic heterozygosity and LD structures for SNPs), we used HAPGEN2 [31–33] to generate genotypes; we calculated SNP MAF and LD structure from 1000 simulated samples. We elected to use the same intersecting set of SNPs for real data and simulation. For HAPGEN2, we eliminated SNPs with MAF<0.002; for 1000 Genomes, we eliminated SNPs for which the call rate (percentage of samples with useful data) was less than 90%. This left n_snp = 11,015,833 SNPs. See S1 Appendix for further details.

We analyzed summary statistics for sixteen phenotypes (S1 DataSourceList.; in what follows, where sample sizes varied by SNP, we quote the median value): (1) major depressive disorder (N_cases = 59,851, N_controls = 113,154) [34]; (2) bipolar disorder (N_cases = 20,352, N_controls = 31,358) [35]; (3) schizophrenia (N_cases = 35,476, N_controls = 46,839) [36]; (4) coronary artery disease (N_cases = 60,801, N_controls = 123,504) [37]; (5) ulcerative colitis (N_cases = 12,366, N_controls = 34,915) and (6) Crohn’s disease (N_cases = 12,194, N_controls = 34,915) [38]; (7) late onset Alzheimer’s disease (LOAD; N_cases = 17,008, N_controls = 37,154) [39] (in the S1 Appendix we present results for a more recent GWAS with N_cases = 71,880 and N_controls = 383,378 [40]); (8) amyotrophic lateral sclerosis (ALS) (N_cases = 12,577, N_controls = 23,475) [41]; (9) number of years of formal education (N = 293,723) [42]; (10) intelligence (N = 262,529) [43, 44]; (11) body mass index (N = 233,554) [45]; (12) height (N = 251,747) [46]; (13) putamen volume (normalized by intracranial volume, N = 11,598) [47]; (14) low- (N = 89,873) and (15) high-density lipoprotein (N = 94,295) [48]; and (16) total cholesterol (N = 94,579) [48]. Most participants were of European ancestry.

For height, we focused on the 2014 GWAS [46], not the more recent 2018 GWAS [49], although we also report below model results for the latter. There are issues pertaining to population structure in the various height GWAS [50, 51], and the 2018 GWAS is a combination of GIANT and UKB GWAS, so some caution is warranted in interpreting results for these data.

For the ALS GWAS data, there is very little signal outside chromosome 9: the data QQ plot essentially tracks the null distribution straight line. The QQ plot for chromosome 9, however, shows a significant departure from the null distribution. Of 471,607 SNPs on chromosome 9 a subset of 273,715 have z-scores, of which 107 are genome-wide significant, compared with 114 across the full genome. Therefore, we restrict ALS analysis to chromosome 9.

A limitation in the current work is that we have not taken account of imputation inaccuracy, where lower MAF SNPs are, through lower LD, less certain. Thus, the effects from lower MAF causal variants will be noisier than for higher MAF variants.

Simulations

We generated genotypes for 10⁵ unrelated simulated samples using HAPGEN2 [33]. For narrow-sense heritability h² equal to 0.1, 0.4, and 0.7, we considered polygenicity π₁ equal to 10⁻⁵, 10⁻⁴, 10⁻³, and 10⁻². For each of these 12 combinations, we randomly selected n_causal = π₁ × n_snp “causal” SNPs and assigned them β-values drawn from the standard normal distribution (i.e., independent of H), with all other SNPs having β = 0. We repeated this ten times, giving ten independent instantiations of random vectors of β’s. Defining Y_G = Gβ, where G is the genotype matrix and β here is the vector of true coefficients over all SNPs, the total phenotype vector is constructed as Y = Y_G+ε, where the residual random vector ε for each instantiation is drawn from a normal distribution such that h² = var(Y_G)/var(Y). For each of the instantiations this implicitly defines the “true” value σ β 2.

The sample simple linear regression slope, β ^, and the Pearson correlation coefficient, r ^, are assumed to be t-distributed. These quantities have the same t-value: t = β ^ / se ( β ^ ) = r ^ / se ( r ^ ) = r ^ N - 2 / 1 - r ^ 2, with corresponding p-value from Student’s t cumulative distribution function (cdf) with N − 2 degrees of freedom: p = 2 × tcdf(−|t|, N − 2) (see S1 Appendix). Since we are not here dealing with covariates, we calculated p from correlation, which is slightly faster than from estimating the regression coefficient. The t-value can be transformed to a z-value, giving the z-score for this p: z = - Φ - 1 ( p / 2 ) × sign ( r ^ ), where Φ is the normal cdf (z and t have the same p-value).

Parameter estimation

We randomly pruned SNPs using the threshold r² > 0.8 to identify “synonymous” SNPs, performing ten such iterations. That is, for each of ten iterations, we randomly selected a SNP (not necessarily the one with largest z-score) to represent each subset of synonymous SNPs. For schizophrenia, for example, pruning resulted in approximately 1.3 million SNPs in each iteration.

The postulated pdf for a SNP’s z-score depends on the SNP’s LD and heterozygosity structure (histogram), H. Given the data–the set of z-scores for available SNPs, as well as their LD and heterozygosity structure—and the H-dependent pdf for z-scores, the objective is to find the model parameters that best predict the distribution of z-scores. We bin the SNPs with respect to a grid of heterozygosity and total LD; for any given H-L bin there will be a range of z-scores whose distribution the model it intended to predict. We find that a 10 × 10 grid of equally spaced bins is adequate for converged results. (Using equally-spaced bins might seem inefficient because of the resulting very uneven distribution of z-scores among grid elements—for example, orders of magnitude more SNPs in grid elements with low total LD compared with high total LD. However, the objective is to model the effects of H and L: using variable grid element sizes so as to maximize balance of SNP counts among grid elements means that the true H- and L-mediated effects of the SNPs in a narrow range of H and L get subsumed with the effects of many more SNPs in a much wider range of H and L—a misspecification of the pdf leading to some inaccuracy.) In lieu of or in addition to total LD (L) binning, one can bin SNPs with respect to their total LD block size (total number of SNPs in LD, ranging from 1 to ≃1,500).

To find the model parameters that best fit the data, for a given H-L bin we binned the selected SNPs z-scores into equally-spaced bins of width dz = 0.0742 (between z_min = −38 and z_max = 38, allowing for p-values near the numerical limit of 10⁻³¹⁶), and from Eq 29 calculated the probability for z-scores to be in each of those z-score bins (the prior probability for “success” in each z-score bin). Then, knowing the actual numbers of z-scores (numbers of “successes”) in each z-score bin, we calculated the multinomial probability, p_m, for this outcome. The optimal model parameter values will be those that maximize the accrual of this probability over all H-L bins. We constructed a cost function by calculating, for a given H-L bin, −ln(p_m) and averaging over prunings, and then accumulating this over all H-L bins. Model parameters minimizing the cost were obtained from Nelder-Mead multidimensional unconstrained nonlinear minimization of the cost function, using the Matlab function fminsearch().

Posterior effect sizes

Model posterior effect sizes, given z (along with N, H, and the model parameters), were calculated using numerical integration over the random variable δ:

Here, since z | δ ∼ N ( δ , σ 0 2 ), the posterior probability of z given δ is simply

P(z) is shorthand for pdf ( z | N , H , π 1 , σ β , σ 0 ), given by Eq 29. P(δ) is calculated by a similar procedure that lead to Eq 29 but ignoring the environmental contributions {E_j}. Specifically, let F_δ = (G₁,…,G_n) denote the vector of products of Fourier transform values. Then, the vector of pdf values for genetic effect bins (indexed by i; numerically, these will be the same as the z-score bins) in the H-L bin, pdf_δ = (f₁,…,f_n) where f i ≡ pdf ( δ i | H ), is

Similarly,

which is used in power calculations.

GWAS replication

A related matter has to do with whether z-scores for SNPs reaching genome-wide significance in a discovery-sample are compatible with the SNPs’ z-scores in a replication-sample, particularly if any of those replication-sample z-scores are far from reaching genome-wide significance, or whether any apparent mismatch signifies some overlooked inconsistency. The model pdf allows one to make a principled statistical assessment in such cases. We present the details for this application, and results applied to studies of bipolar disorder, in the S1 Appendix (pp. S7-S11).

GWAS power

Chip heritability, h S N P 2, is the proportion of phenotypic variance that in principle can be captured additively by the n_snp SNPs under study [17]. It is of interest to estimate the proportion of h S N P 2 that can be explained by SNPs reaching genome-wide significance, p≤5 × 10⁻⁸ (i.e., for which |z|>z_t = 5.45), at a given sample size [52, 53]. In Eq 1, for SNP i with genotype vector g_i over N samples, let y g i ≡ g i β i. If the SNP’s heterozygosity is H_i, then var ( y g i )=β i 2 H i. If we knew the full set {β_i} of true β-values, then, for z-scores from a particular sample size N, the proportion of SNP heritability captured by genome-wide significant SNPs, A(N), would be given by

Now, from Eq 15, δ i = N ∑ j H j r i j β j. If SNP i is causal and sufficiently isolated so that it is not in LD with other causal SNPs, then δ i = N H i β i, and var ( y g i ) = δ i 2 / N. When all causal SNPs are similarly isolated, Eq 34 becomes

Of course, the true β_i are not known and some causal SNPs will likely be in LD with others. Furthermore, due to LD with causal SNPs, many SNPs will have a nonzero (latent or unobserved) effect size, δ. Nevertheless, we can formulate an approximation to A(N) which, assuming the pdf for z-scores (Eq 29) is reasonable, will be inaccurate to the degree that the average LD structure of genome-wide significant SNPs differs from the overall average LD structure. As before (see the subsection “Model PDF: Convolution”), consider a fixed set of n equally-spaced nominal z-scores covering a wide range of possible values (changing from the summations in Eq 35 to the uniform summation spacing Δz now requires bringing the probability density into the summations). For each z from the fixed set (and, as before, employing data reduction by averaging so that H and L denote values for the 10 × 10 grid), use E(δ²|z, N, H, L) given in Eq 33 to define

(emphasizing dependence on N, H, and L). Then, for any N, A(N) can be estimated by

where ∑_H,L denotes sum over the H-L grid elements. The ratio in Eq 37 should be accurate if the average effects of LD in the numerator and denominator cancel—which will always be true as the ratio approaches 1 for large N. Plotting A(N) gives an indication of the power of future GWAS to capture chip heritability.

Quantile-quantile plots and genomic control

One of the advantages of quantile-quantile (QQ) plots is that on a logarithmic scale they emphasize behavior in the tails of a distribution, and provide a valuable visual aid in assessing the independent effects of polygenicity, strength of association, and variance distortion—the roles played by the three model parameters–as well as showing how well a model fits data. QQ plots for the model were constructed using Eq 29, replacing the normal pdf with the normal cdf, and replacing z with an equally-spaced vector z → n o m of length 10,000 covering a wide range of nominal |z| values (0 through 38). SNPs were divided into a 10 × 10 grid of H × L bins, and the cdf vector (with elements corresponding to the z-values in z → n o m) accumulated for each such bin (using mean values of H and L for SNPs in a given bin).

For a given set of samples and SNPs, the genomic control factor, λ, for the z-scores is defined as the median z² divided by the median for the null distribution, 0.455 [19]. This can also be calculated from the QQ plot. In the plots we present here, the abscissa gives the -log₁₀ of the proportion, q, of SNPs whose z-scores exceed the two-tailed significance threshold p, transformed in the ordinate as -log_10(p). The median is at q_med = 0.5, or −log₁₀(q_med) ≃ 0.3; the corresponding empirical and model p-value thresholds (p_med) for the z-scores—and equivalently for the z-scores-squared—can be read off from the plots. The genomic inflation factor is then given by

Note that the values of λ reported here are for pruned SNP sets; these values will be lower than for the total GWAS SNP sets.

Knowing the total number, n_tot, of p-values involved in a QQ plot (number of GWAS z-scores from pruned SNPs), any point (q, p) (log-transformed) on the plot gives the number, n_p = qn_tot, of p-values that are as extreme as or more extreme than the chosen p-value. This can be thought of as n_p “successes” out of n_tot independent trials (thus ignoring LD) from a binomial distribution with prior probability q. To approximate the effects of LD, we estimate the number of independent SNPs as n_tot/f where f ≃ 10. The 95% binomial confidence interval for q is calculated as the exact Clopper-Pearson 95% interval [54], which is similar to the normal approximation interval, q ± 1 . 96 q ( 1 - q ) / n t o t / f.

Number of causal SNPs

The estimated number of causal SNPs is given by the polygenicity, π₁, times the total number of SNPs, n_snp: n_causal = π₁ n_snp. n_snp is given by the total number of SNPs that went into building the heterozygosity/LD structure, H in Eq 29, i.e., the approximately 11 million SNPs selected from the 1000 Genomes Phase 3 reference panel, not the number of typed SNPs in the particular GWAS. The parameters estimated are to be seen in the context of the reference panel, which we assume contains all common causal variants. Stable quantities (i.e., fairly independent of the reference panel size. e.g., using the full panel or ignoring every second SNP), are the estimated effect size variance and number of causal variants—which we demonstrate below—and hence the heritability. Thus, the polygenicity will scale inversely with the reference panel size. A reference panel with a substantially larger number of samples would allow for inclusion of more SNPs (non-zero MAF), and thus the actual polygenicity estimated would change slightly.

Narrow-sense chip heritability

Since we are treating the β coefficients as fixed effects in the simple linear regression GWAS formalism, with the phenotype vector standardized with mean zero and unit variance, from Eq 1 the proportion of phenotypic variance explained by a particular causal SNP whose reference panel genotype vector is g, q² = var(y;g), is given by q² = β² H. The proportion of phenotypic variance explained additively by all causal SNPs is, by definition, the narrow sense chip heritability, h². Since E(β²)=σ β 2 and n_causal = π₁ n_snp, and taking the mean heterozygosity over causal SNPs to be approximately equal to the mean over all SNPs, H ¯, the chip heritability can be estimated as

Mean heterozygosity from the ≃11 million SNPs is H ¯ = 0 . 2165.

For all-or-none traits like disease status, the estimated h² from Eq 38 for an ascertained case-control study is on the observed scale and is a function of the prevalence in the adult population, K, and the proportion of cases in the study, P. The heritability on the underlying continuous liability scale [55], h l 2, is obtained by adjusting for ascertainment (multiplying by K(1 − K)/(P(1 − P)), the ratio of phenotypic variances in the population and in the study) and rescaling based on prevalence [6, 56]:

where a is the height of the standard normal pdf at the truncation point z_K defined such that the area under the curve in the region to the right of z_K is K.

Confidence intervals

Confidence intervals for parameters were estimated using the inverse of the observed Fisher information matrix (FIM). The full FIM was estimated for all three parameters used in the model. For the derived quantity h², which depends on all parameters, the covariances among the parameters, given by the off-diagonal elements of the inverse of the FIM, were incorporated. Numerical values are in S1 Appendix (p. S12).

Results

Simulations

Table 1 shows the simulation results, comparing true and estimated values for the model parameters, heritability, and the number of causal SNPs, for twelve scenarios where π₁ and σ β 2 both range over three orders of magnitude, encompassing the range of values for the phenotypes; in S1 Appendix (p. S15) are QQ plots for a randomly chosen (out of 10) β-vector and phenotype instantiation for each of the twelve (π₁, h²) scenarios. Most of the π ^ 1 estimates are in very good agreement with the true values, though for the extreme scenario of high heritability and low polygenicity it is overestimated by factors of two-to-three. The numbers of estimated causal SNPs (out of ≃11 million) are in correspondingly good agreement with the true values, ranging in increasing powers of 10 from 110 through 110,158. The estimated discoverabilities (σ ^ β 2) are also in good agreement with the true values. In most cases, σ ^ 0 2 is close to 1, indicating little or no global inflation, though it is elevated for high heritability with high polygenicity, suggesting it is capturing some ubiquitous effects.

In the S1 Appendix (pp. S5-S7) we examine the issue of model misspecification. Specifically, we assign causal effects β drawn from a Gaussian whose variance is not simply a constant but depends on heterozygosity, such that rarer causal SNPs will tend to have larger effects [15]. The results—see S1 Appendix (p. S5)–show that the model still makes reasonable estimates of the underlying genetic architecture. Additionally, we tested the scenario where true causal effects are distributed with respect to two Gaussians [14], a situation that allows for a small number of the causal SNPs to have quite large effects—see S1 Appendix (p. S6). We find that heritabilities are still reasonably estimated using our model. In all these scenarios the overall data QQ plots were accurately reproduced by the model. As a counter example, we simulated summary statistics where the prior probability of a reference SNP being causal decreased linearly with total LD (see S1 Appendix (p. S7)). In this case, our single Gaussian fit (which assumes no LD dependence on the prior probability of a reference SNP being causal) did not produce model QQ plots that accurately tracked the data QQ plots (see S1 Appendix (p. S19)). The model parameters and heritabilities were also poor. But this scenario is highly artificial; in contrast, in situations where the data QQ plots were accurately reproduced by the model, the estimated model parameters and heritability were plausible.

Phenotypes

Figs 1 and 2 show QQ plots for the pruned z-scores for eight qualitative and eight quantitative phenotypes, along with model estimates (S1 Appendix (pp. S20-S38) show a 4 × 4 grid breakdown with respect to heterozygosity × total-LD of QQ plots for all phenotypes studied here; the 4 × 4 grid is a subset of the 10 × 10 grid used in the calculations). In all cases, the model fit (yellow) closely tracks the data (dark blue). For the sixteen phenotypes, estimates for the model polygenicity parameter (fraction of reference panel, with ≃11 million SNPs, estimated to have non-null effects) range over two orders of magnitude, from π₁ ≃ 2 × 10⁻⁵ to π₁ ≃ 4 × 10⁻³. The estimated SNP discoverability parameter (variance of β, or expected β², for causal variants) also ranges over two orders of magnitude from σ β 2 ≃ 7 × 10 - 6 to σ β 2 ≃ 2 × 10 - 3 (in units where the variance of the phenotype is normalized to 1).

We find that schizophrenia and bipolar disorder appear to be similarly highly polygenic, with model polygenicities ≃ 2.84 × 10⁻³ and ≃ 2.70 × 10⁻³, respectively. The model polygenicity of major depressive disorder, however, is 40% higher, π₁ ≃ 4 × 10⁻³—the highest value among the sixteen phenotypes. In contrast, the model polygenicities of late onset Alzheimer’s disease and Crohn’s disease are almost thirty times smaller than that of schizophrenia.

In S1 Appendix (p. S17) we show results for Alzheimer’s disease exclusively for chromosome 19 (which contains APOE), and for all autosomal chromosomes excluding chromosome 19. We also show results with the same chromosomal breakdown for a recent GWAS involving 455,258 samples that included 24,087 clinically diagnosed LOAD cases and 47,793 AD-by-proxy cases (individuals who were not clinically diagnosed with LOAD but for whom at least one parent had LOAD) [65]. These GWAS give consistent estimates of polygenicity: π₁ ∼ 1 × 10⁻⁴ excluding chromosome 19, and π₁ ∼ 6 × 10⁻⁵ for chromosome 19 exclusively.

Of the quantitative traits, educational attainment has the highest model polygenicity, π₁ = 3.2 × 10⁻³, similar to intelligence, π₁ = 2.2 × 10⁻³. Approximately two orders of magnitude lower in polygenicity are the endophenotypes putamen volume and low- and high-density lipoproteins.

The model effective SNP discoverability for schizophrenia is σ ^ β 2 = 5 . 51 × 10 - 5, similar to that for bipolar disorder. Major depressive disorder, which has the highest polygenicity, has the lowest SNP discoverability, approximately one-eighth that of schizophrenia; it is this low value, combined with high polygenicity that leads to the weak signal in Fig 1 (A) even though the sample size is relatively large. In contrast, SNP discoverability for Alzheimer’s disease is almost four times that of schizophrenia. The inflammatory bowel diseases, however, have much higher SNP discoverabilities, 16 and 31 times that of schizophrenia respectively for ulcerative colitis and Crohn’s disease—the latter having the second highest value of the sixteen phenotypes: σ ^ β 2 = 1 . 7 × 10 - 3.

Additionally, for Alzheimer’s disease we show in the S1 Appendix (p. S17) that the discoverability is two orders of magnitude greater for chromosome 19 than for the remainder of the autosome. Note that since two-thirds of the 2018 “cases” are AD-by-proxy, the discoverabilities for the 2018 data are, as expected, reduced relative to the values for the 2013 data (approximately 3.5 times smaller).

The narrow sense SNP heritability from the ascertained case-control schizophrenia GWAS is estimated as h² = 0.37. Taking adult population prevalence of schizophrenia to be K = 0.01 [66, 67] (but see also [68], for K = 0.005), and given that there are 51,900 cases and 71,675 controls in the study, so that the proportion of cases in the study is P = 0.42, the heritability on the liability scale for schizophrenia from Eq 39 is h ^ l 2 = 0 . 21. For bipolar disorder, with K = 0.005 [58], 20,352 cases and 31,358 controls, h ^ l 2=0.16. Major depressive disorder appears to have a much lower model-estimated SNP heritability than schizophrenia: h ^ l 2 = 0.07. The model estimate of SNP heritability for height is 17%, lower than the oft-reported value ≃50% (see Discussion). However, despite the huge differences in sample size, we find the same value, 17%, for the 2010 GWAS (N = 133,735 [69]), and 19% for the 2018 GWAS (N = 707,868 [46, 49])—see Table 2.

Fig 3 shows the sample size required so that a given proportion of chip heritability is captured by genome-wide significant SNPs for the phenotypes (assuming equal numbers of cases and controls for the qualitative phenotypes: N_eff = 4/(1/N_cases + 1/N_controls), so that when N_cases = N_controls, N_eff = N_cases + N_controls = N, the total sample size, allowing for a straightforward comparison with quantitative traits). At current sample sizes, only 4% of narrow-sense chip heritability is captured for schizophrenia and only 1% for bipolar disorder; using current methodologies, a sample size of N_eff ∼ 1 million would be required to capture the preponderance of SNP heritability for these phenotypes. Major depressive disorder GWAS currently is greatly under-powered, as shown in Fig 3(A). For education, we predict that 3.5% of phenotypic variance would be explained at N = 1.1 million, in good agreement with the value found from direct computation of 3.2% [70]. For other phenotypes, the proportions of total SNP heritability captured at the available sample sizes are given in Fig 3.

The sample size for ALS was quite low, and we restricted the analysis to chromosome 9, which had most of the genome-wide significant typed SNPs; we estimate that there are ≃7 causal SNPs with high discoverability on chromosome 9 [71, 72], with very high discoverability, σ β 2 ≃ 0.003. In contrast, for AD restricted to chromosome 19, there were an estimated 14 causal SNPs with discoverability σ β 2 ≃ 0.02 (see the S1 Appendix (p. S17)).

In this study, we assume that population stratification in the raw data has been corrected for in the publicly-available summary statistics. However, given that some of the sample sizes are extremely large, we allow for the possibility of residual cryptic relatedness. This would result in a scaling of the z-scores, Eq 9 [19]. Thus, to test the modeling of inflation due to cryptic relatedness, we scaled the simulation z-scores as described earlier (z = σ₀ z_u with σ₀ > 1, where z_u are the original z-scores, i.e., not artificially inflated) and reran the model. E.g., for education and schizophrenia we inflated the z-scores by a factor of 1.2. For schizophrenia we found σ 0 2 = 1 . 366, which is almost exactly as predicted (1.14 × 1.2 = 1.368), while the polygenicity and discoverability parameters are essentially unchanged: π₁ = 2.81 × 10⁻³, and σ β 2 = 5 . 56 × 10 - 5. For education we found σ 0 2 = 1 . 206, which again is almost exactly as predicted (1.0 × 1.2 = 1.2), while the polygenicity and discoverability parameters are again essentially unchanged: π₁ = 3.19 × 10⁻³, and σ β 2 = 1 . 57 × 10 - 5.

A comparison of our results with those of [14] and [15] is in the S1 Appendix (p. S13). Critical methodological differences with model M2 in [14] are that we use a full reference panel of 11 million SNPs from 1000 Genomes Phase 3, we allow for the possibility of inflation in the data, and we provide an exact solution, based on Fourier Transforms, for the z-score pdf arising from the posited distribution of causal effects, resulting in better fits of the model and the data QQ plots—as can be seen by comparing our QQ plots with those reported in S1 Appendix (p. S13). Although our estimated number of causal are often within a factor of two of those from the nominally equivalent model M2 of Zhang et al, there is no clear pattern to the mismatch.

GWAS replication

In the S1 Appendix (pp. S7-S11) we provide an extensive example of testing the compatibility of summary statistics from two large bipolar disorder GWASs. Because z-scores are so noisy, it is possible for a typed SNP with a highly significant p-values in one GWAS to completely fail to reach significance in a subsequent GWAS, and for these outcomes to be statistically consistent. SNP heterozygosity and total LD, as well as sample sizes, are relevant in making such assessments.

Dependence on reference panel

Given a liberal MAF threshold of 0.002, our reference panel should contain the vast majority of common SNPs for European ancestry. However, it does not include other structural variants (such as small insertions/deletions, or haplotype blocks) which may also be causal for phenotypes. To validate our parameter estimates for an incomplete reference, we reran our model on real phenotypes using the culled reference where we exclude every other SNP. The result is that all estimated parameters are as before except that π ^ 1 doubles, leaving the estimatde number of causal SNPs and heritability as before. For example, for schizophrenia we get π₁ = 5.3 × 10⁻³ and σ β 2 = 5 . 8 × 10 - 5 for the reduced reference panel, versus π₁ = 2.8 × 10⁻³ and σ β 2 = 5 . 5 × 10 - 5 for the full panel, with heritability remaining essentially the same (37% on the observed scale).

Discussion

Here we present a unified method based on GWAS summary statistics, incorporating detailed LD structure from an underlying reference panel of SNPs with MAF>0.002, for estimating: phenotypic polygenicity, π₁, expressed as the fraction of the reference panel SNPs that have a non-null true β value, i.e., are “causal”; and SNP discoverability or mean strength of association (the variance of the underlying causal effects), σ β 2. In addition the model can be used to estimate residual inflation of the association statistics due to variance distortion induced by cryptic relatedness, σ 0 2. The model assumes that there is very little, if any, inflation in the GWAS summary statistics due to population stratification (bias shift in z-scores due to ethnic variation).

We apply the model to sixteen diverse phenotypes, eight qualitative and eight quantitative. From the estimated model parameters we also estimate the number of causal common-SNPs in the underlying reference panel, n_causal, and the narrow-sense common-SNP heritability, h² (for qualitative phenotypes, we re-express this as the proportion of population variance in disease liability, h l 2, under a liability threshold model, adjusted for ascertainment); in the event rare SNPs (i.e., not in the reference panel) are causal, h² will be an underestimate of the true SNP heritability. In addition, we estimate the proportion of SNP heritability captured by genome-wide significant SNPs at current sample sizes, and predict future sample sizes needed to explain the preponderance of SNP heritability.

We find that schizophrenia is highly polygenic, with π₁ = 2.8 × 10⁻³. This leads to an estimate of n_causal ≃ 31, 000, which is in reasonable agreement with a recent estimate that the number of causal SNPs is >20,000 [73]. The SNP associations, however, are characterized by a narrow distribution, σ β 2 = 6 . 27 × 10 - 5, indicating that most associations are of weak effect, i.e., have low discoverability. Bipolar disorder has similar parameters. The smaller sample size for bipolar disorder has led to fewer SNP discoveries compared with schizophrenia. However, from Fig 3, sample sizes for bipolar disorder are approaching a range where rapid increase in discoveries becomes possible. For educational attainment [42, 74, 75], the polygenicity is somewhat greater, π₁ = 3.2 × 10⁻³, leading to an estimate of n_causal ≃35, 000, half a recent estimate, ≃ 70, 000, for the number of loci contributing to heritability [74]. The variance of the distribution for causal effect sizes is a quarter that of schizophrenia, indicating lower discoverability. Intelligence, a related phenotype [43, 76], has a larger discoverability than education while having lower polygenicity (∼ 10, 000 fewer causal SNPs).

In marked contrast are the lipoproteins and putamen volume which have very low polygenicity: π₁ < 5 × 10⁻⁵, so that only 250 to 550 SNPs (out of ≃11 million) are estimated to be causal. However, causal SNPs for putamen volume and HDL appear to be characterized by relatively high discoverability, respectively 17-times and 23-times larger than for schizophrenia (see S1 Appendix (p. S9) for additional notes on HDL, and [77] for a relevant comparison with our work).

The QQ plots (which are sample size dependent) reflect these differences in genetic architecture. For example, the early departure of the schizophrenia QQ plot from the null line indicates its high polygenicity, while the steep rise for putamen volume after its departure corresponds to its high SNP discoverability.

For Alzheimer’s disease, our estimate of the liability-scale SNP heritability for the full 2013 dataset [39] is 15% for prevalence of 14% for those aged 71 older, half from APOE, while the recent “M2” and “M3” models of Zhang et al [14] gave values of 7% and 10% respectively–see S1 Appendix (p. S13). A recent report from two methods, LD Score Regression (LDSC) and SumHer [77], estimated SNP heritability of 3% for LDSC and 12% for SumHer (assuming prevalence of 7.5%). A raw genotype-based analysis (GCTA), including genes that contain rare variants that affect risk for AD, reported SNP heritability of 53% [7, 78]; an earlier related study that did not include rare variants and had only a quarter of the common variants estimated SNP heritability of 33% for prevalence of 13% [79]. GCTA calculations of heritability are within the domain of the so-called infinitesimal model where all markers are assumed to be causal. Our model suggests, however, that phenotypes are characterized by polygenicities less than 5 × 10⁻³; for AD the polygenicity is ≃ 10⁻⁴. Nevertheless, the GCTA approach yields a heritability estimate closer to the twin-based (broad sense) value, estimated to be in the range 60-80% [80]. The methodology appears to be robust to many assumptions about the distribution of effect sizes [81, 82]; the SNP heritability estimate is unbiased, though it has larger standard error than methods that allow for only sparse causal effects [69, 83]. For the 2013 data analyzed here [39], a summary-statistics-based method applied to a subset of 54,162 of the 74,046 samples gave SNP heritability of almost 7% on the observed scale [12, 84]; our estimate is 12% on the observed scale—see S1 Appendix (p. S17).

Onset and clinical progression of sporadic Alzheimer’s disease is strongly age-related [85, 86], with prevalence in differential age groups increasing at least up through the early 90s [62]. Thus, it would be more accurate to assess heritability (and its components, polygenicity and discoverability) with respect to, say, five-year age groups beginning with age 65 years, and using a consistent control group of nonagenarians and centenarians. By the same token, comparisons among current and past AD GWAS are complicated because of potential differences in the age distributions of the respective case and the control cohorts. Additionally, the degree to which rare variants are included will affect heritability estimates. The summary-statistic-based estimates of polygenicity that we report here are, however, likely to be robust for common SNPs: π₁ ≃ 1.1 × 10⁻⁴, with only a few causal SNPs on chromosome 19.

Our point estimate for the liability-scale SNP heritability of schizophrenia is h l 2 = 0 . 21 (assuming a population risk of 0.01), and that 4% of this (i.e., 1% of overall disease liability) is explainable based on common SNPs reaching genome-wide significance at the current sample size. This h l 2 estimate is in reasonable agreement with a recent result, h l 2 = 0 . 27 [73, 87], also calculated from the PGC2 data set but using raw genotype data for 472,178 markers for a subset of 22,177 schizophrenia cases and 27,629 controls of European ancestry; and with an earlier result of h l 2 = 0 . 23 from PGC1 raw genotype data for 915,354 markers for 9,087 schizophrenia cases and 12,171 controls [7, 88]. The recent “M2” (single non-null Gaussian) model estimate is h l 2 = 0 . 29 [14] (see S1 Appendix (p. S13)). No QQ plot was available for the M2 model fit to schizophrenia data, but such plots (truncated on the y-axis at −log₁₀(p) = 10) for many other phenotypes were reported [14]. We note that for multiple phenotypes (height, LDL cholesterol, total cholesterol, years of schooling, Crohn’s disease, coronary artery disease, and ulcerative colitis) our single causal Gaussian model appears to provide a better fit to the data than M2: many of the M2 plots show a very early and often dramatic deviation between prediction and data, as compared with our model QQ plots which are also built from a single causal Gaussian, suggesting an upward bias in polygenicity and/or variance of effect sizes, and hence heritability as measured by the M2 model for these phenotypes. The LDSC liability-scale (1% prevalence) SNP heritability for schizophrenia has been reported as h l 2 = 0 . 555 [12] and more recently as 0.19 [77], the latter in very good agreement with our estimate; on the observed scale it has been reported as 45% [12, 84], in contrast to our corresponding value of 37%. Our estimate of 1% of overall variation on the liability scale for schizophrenia explainable by genome-wide significant loci compares reasonably with the proportion of variance on the liability scale explained by Risk Profile Scores (RPS) reported as 1.1% using the “MGS” sample as target (the median for all 40 leave-one-out target samples analyzed is 1.19%—see Extended Data Figure 5 and Supplementary Tables 5 and 6 in [36]; this was incorrectly reported as 3.4% in the main paper). These results show that current sample sizes need to increase substantially in order for RPSs to have predictive utility, as the vast majority of associated SNPs remain undiscovered. Our power estimates indicate that ∼500,000 cases and an equal number of controls would be needed to identify these SNPs (note that there is a total of approximately 3 million cases in the US alone).

A subtle but important issue is downward bias of large-sample maximum-likelihood estimates of SNP heritability, due to over-ascertainment of cases in case-control studies [87]; it has been examined in the context of restricted maximum likelihood (REML) in GCTA, which assumes a polygenicity of 1, i.e., every SNP is causal. For schizophrenia, this has been assessed in the context of BOLT-REML, which assumes a mixture distribution of small (‘spike’) and large (‘slab’) effects [73]: from 22,177 cases and 27,629 controls, the observed-scale heritability is reported as h² = 0.415, equivalent to h l 2 = 0 . 23 on the liability scale, assuming 1% disease prevalence. However, using “phenotype correlation-genetic correlation” (PCGC) regression, a moments-based approach requiring raw-genotype data which produces unbiased estimates for case-control studies of disease traits [87], the unbiased liability-scale heritability is reported as h g 2 = 0 . 27, indicating that the likelihood-maximization estimate is biased down by 15% of the unbiased value (the degree of underestimation decreases for smaller sample sizes). Our estimate for the liability-scale heritability of schizophrenia, from a larger sample than in [73], is h l 2 = 0 . 21. This at least would be consistent with downward bias operating in point-normal causal distributions, in a manner similar to that in GCTA and BOLT-REML. This would then translate into either an underestimate of the number of causal SNPs, or more likely an underestimate of the variance of the distribution of causal effects.

For educational attainment, we estimate SNP heritability h² = 0.12, in good agreement with the estimate of 11.5% given in [42]. As with schizophrenia, this is substantially less than the estimate of heritability from twin and family studies of ≃40% of the variance in educational attainment explained by genetic factors [74, 89].

For putamen volume, we estimate the SNP heritability h² = 0.11, in reasonable agreement with an earlier estimate of 0.1 for the same overall data set [4, 47]. For LDL and HDL, we estimate h² = 0.06 and h² = 0.07 respectively, in good agreement with the LDSC estimates h² = 0.08 and h² = 0.07 [77], and the M2 model of [14]—see S1 Appendix (p. S13).

For height (N = 251,747 [46]) we find that its model polygenicity is π₁ = 5.66 × 10⁻⁴, a quarter that of intelligence, while its discoverability is five times that of intelligence, leading to a SNP heritability of 17%. The number of causal SNPs (out of a total of approximately 11 million) is approximately 6k; although this is about one twentieth the estimate reported in [90], it remains large and allows for height to be interpreted as “omnigenic”. For the 2010 GWAS (N = 133,735 [69]) and 2018 GWAS (N = 707,868 [49]), we estimate SNP heritability of 17% and 19% respectively (see Table 2, and S1 Appendix (p. S17)). These heritabilities are in considerable disagreement with the SNP heritability estimate of ≃50% [46] (average of estimates from five cohorts ranging in size from N = 1,145 to N = 5,668, with ≃1 million SNPs). For the 2010 GWAS, the M2 model [14] gives h² = 0.30 (see S1 Appendix (p. S13)); the upward deviation of the model QQ plot in [14] suggests that this value might be inflated. For the 2014 GWAS, the M3 model estimate is h² = 33% [14]; the Regression with Summary Statistics (RSS) model estimate is h² = 52% (with ≃11, 000 causal SNPs) [91], which, not taking any inflation into account, is definitely a model overestimate; and in [77] the LDSC estimate is reported as h² = 20% while the SumHer estimate is h² = 46% (in general across traits, the SumHer heritability estimates tend to be two-to-five times larger than the LDSC estimates). The M2, M3, and RSS models use a reference panel of ≃1 million common SNPs, in contrast with the ≃11 million SNPs used in our analysis. Also, it should be noted that the M2, M3, and RSS model estimates did not take the possibility of inflation into account. For the 2014 height GWAS, that inflation is reported as the LDSC intercept is 2.09 in [77], indicating considerable inflation; for the 2018 dataset we find σ 0 2 = 2.5, while the LD score regression intercept is 2.1116 (se 0.0458). Given the various estimates of inflation and the controversy over population structure in the height data [50, 51], it is not clear what results are definitely incorrect.

Our power analysis for height (2014) shows that 37% of the narrow-sense heritability arising from common SNPs is explained by genome-wide significant SNPs (p ≤ 5 × 10⁻⁸), i.e., 6.3% of total phenotypic variance, which is substantially less than the 16% direct estimate from significant SNPs [46]. It is not clear why these large discrepancies exist. One relevant factor, however, is that we estimate a considerable confounding (σ 0 2 = 1.66) in the height 2014 dataset. Our h² estimates are adjusted for the potential confounding measured by σ 0 2, and thus they represent what is likely a lower bound of the actual SNP-heritability, leading to a more conservative estimate than what has previously been reported. We note that after adjustment, our h² estimates are consistent across all three datasets (height 2010, 2014 and 2018), which otherwise would range by more than 2.5-fold. Another factor might be the relative dearth of typed SNPs with low heterozygosity and low total LD (see top left segment in S1 Appendix (p. S29), n = 780): there might be many causal variants with weak effect that are only weakly tagged. Nevertheless, given the discrepancies noted above, caution is warranted in interpreting our model results for height.

Conclusion

The common-SNP causal effects model we have presented is based on GWAS summary statistics and detailed LD structure of an underlying reference panel, and assumes a Gaussian distribution of effect sizes at a fraction of SNPs randomly distributed across the autosomal genome. While not incorporating the effects of rare SNPs, we have shown that it captures the broad genetic architecture of diverse complex traits, where polygenicities and the variance of the effect sizes range over orders of magnitude.

The current model (essentially Eq 4) and its implementation (essentially Eq 29) are basic elements for building a more refined model of SNP effects using summary statistics. Higher accuracy in characterizing causal alleles in turn will enable greater power for SNP discovery and phenotypic prediction.

Supporting information

S1 Appendix [pdf]
Appendix.

S1 DataSourceList [xlsx]
Data source list.