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, (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) (1) 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 (2) where T denotes transpose and 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, : (3)

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] (4) where is the Dirac delta function, so that considering all SNPs, the net variance is . 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, , which is assumed to be independent of δ [18]: (5) with (6) so that (7) where (8)

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 (9) where σ₀ ≥ 1 characterizes the inflation. Thus, from Eq 7, in the presence of inflation in the form of variance distortion (10) where , so that and . In the presence of variance distortion one is dealing with inflated random variables , 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 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., . So a more correct δ is (11)

Taylor expanding in q and then taking the variance gives (12)

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, , 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 (13) (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; 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 is simply cov(g_i, g_j), the covariance between genotypes for SNPs i and j. Since correlation is covariance normalized by the variances, in Eq 13 can be written as (14) 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: (15)

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, [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 in Eq 10 will be given by (16)

For this idealized case, the marginal distribution, or pdf, of z-scores for a set of such associated SNPs is (17) where ϕ(⋅;μ, σ²) is the normal distribution with mean μ and variance σ², and 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 (18) 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 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 and 1 where (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 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 , 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 () 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 ( 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 (19) and the corresponding Gaussian is (20)

For a SNP with LD and heterozygosity structure , 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: (21) where K_max is bounded above by n. Note again that 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 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 , 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, (22)

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 . 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), (23)

The Fourier transform (with variable k—see below) of the first term on the right hand side is (24) while that of the second term is simply (1 − π₁). Additionally, the environmental term is (ignoring LD-induced correlation, as noted earlier), and its Fourier transform is . 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 (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; 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 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 (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 r² bin whose average heterozygosity is H_w. Thus, one can re-express z-scores for an H-L bin as (25) 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 , 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 (26)

(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 (27)

Recall that 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 denote the set of model parameters. The Fourier transform of the environmental contribution, denoted E_j ≡ E(k_j), is (28)

Let F_z = (G₁ E₁,…,G_nE_n) denote the vector of products of Fourier transform values, and let 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 , is (29)

Thus, the i^th element 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 .

The sample simple linear regression slope, , and the Pearson correlation coefficient, , are assumed to be t-distributed. These quantities have the same t-value: , 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: , 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), . Given the data–the set of z-scores for available SNPs, as well as their LD and heterozygosity structure—and the -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, , and the model parameters), were calculated using numerical integration over the random variable δ: (30)

Here, since , the posterior probability of z given δ is simply (31) P(z) is shorthand for , 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 , is (32)

Similarly, (33) 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, , 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 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 . If the SNP’s heterozygosity is H_i, then . 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 (34)

Now, from Eq 15, . If SNP i is causal and sufficiently isolated so that it is not in LD with other causal SNPs, then , and = . When all causal SNPs are similarly isolated, Eq 34 becomes (35)

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 (36) (emphasizing dependence on N, H, and L). Then, for any N, A(N) can be estimated by (37) 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 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 ) 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, .

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, 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(β²)= and n_causal = π₁ n_snp, and taking the mean heterozygosity over causal SNPs to be approximately equal to the mean over all SNPs, , the chip heritability can be estimated as (38)

Mean heterozygosity from the ≃11 million SNPs is .

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], , 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]: (39) 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).

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 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 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 () are also in good agreement with the true values. In most cases, 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Simulation results: Comparison of mean (std) true and estimated (^) model parameters and derived quantities.

Results for each line, for specified heritability h² and fraction π₁ of causal SNPs, are from 10 independent instantiations with random selection of the n_causal causal SNPs that are assigned a β-value from the standard normal distribution. Defining Y_g = Gβ, where G is the genotype matrix, 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 var(Y) = var(Y_g)/h² for predefined h². For each of the instantiations, i, this implicitly defines the true value , and is their mean. An example QQ plot for each line entry is shown in in S1 Appendix (p. S15).

https://doi.org/10.1371/journal.pgen.1008612.t001

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 to (in units where the variance of the phenotype is normalized to 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. QQ plots of (pruned) z-scores for qualitative phenotypes (dark blue, 95% confidence interval in light blue) with model prediction (yellow): (A) major depressive disorder; (B) bipolar disorder; (C) schizophrenia; (D) coronary artery disease (CAD); (E) ulcerative colitis (UC); (F) Crohn’s disease (CD); (G) late onset Alzheimer’s disease (AD), excluding APOE (see also S1 Appendix (p. S17)); and (H) amyotrophic lateral sclerosis (ALS), restricted to chromosome 9 (see also S1 Appendix (p. S18)).

The dashed line is the expected QQ plot under null (no SNPs associated with the phenotype). p is a nominal p-value for z-scores, and q is the proportion of z-scores with p-values exceeding that threshold. λ is the overall nominal genomic control factor for the pruned data (which is accurately predicted by the model in all cases). The three estimated model parameters are: polygenicity, ; discoverability, (corrected for inflation); and SNP association χ²-statistic inflation factor, . is the estimated narrow-sense chip heritability, re-expressed as on the liability scale for these case-control conditions assuming a prevalence of: MDD 7.1% [57], BIP 0.5% [58], SCZ 1% [59], CAD 3% [60], UC 0.1% [61], CD 0.1% [61], AD 14% (for people aged 71 and older in the USA [62, 63]), and ALS 5 × 10⁻⁵ [64]. The estimated number of causal SNPs is given by where n_snp = 11, 015, 833 is the total number of SNPs, whose LD structure and MAF underlie the model; the GWAS z-scores are for subsets of these SNPs. N_eff is the effective case-control sample size–see text. Reading the plots: on the vertical axis, choose a p-value threshold (more extreme values are further from the origin), then the horizontal axis gives the proportion of SNPs exceeding that threshold (higher proportions are closer to the origin). Numerical values for the model parameters are also given in Table 2. See also S1 Appendix (pp. S20-S28).

https://doi.org/10.1371/journal.pgen.1008612.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. this area for notes QQ plots of (pruned) z-scores and model fits for quantitative phenotypes: (A) educational attainment; (B) intelligence; (C) body mass index (BMI); (D) height; (E) putamen volume; (F) low-density lipoprotein (LDL); (G) high-density lipoprotein (HDL); and (H) total cholesterol (TC).

N is the sample size. See Fig 1 for further description. Numerical values for the model parameters are also given in Table 2. See also S1 Appendix (pp. S29-S38).

https://doi.org/10.1371/journal.pgen.1008612.g002

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 , 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: .

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 . For bipolar disorder, with K = 0.005 [58], 20,352 cases and 31,358 controls, . Major depressive disorder appears to have a much lower model-estimated SNP heritability than schizophrenia: = 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of model results for phenotypes shown in Figs 1 and 2.

The subscript in indicates that for the qualitative phenotypes (the first eight) the reported SNP heritability is on the liability scale. MDD: Major Depressive Disorder; CAD: coronary artery disease; AD: Alzheimer’s Disease (excluding APOE locus; *for the full autosomal reference panel, i.e., including APOE, for AD—see S1 Appendix (p. S17)); BMI: body mass index; ^†ALS: amyotrophic lateral sclerosis, restricted to chromosome 9; LDL: low-density lipoproteins; HDL: high-density lipoproteins. ^$In addition to the 2014 height GWAS (N = 251,747 [46]), we include here model results for the 2010 (N = 133,735 [69]) and 2018 (N = 707,868 [49]) height GWAS; there is remarkable consistency for the 2010 and 2014 GWAS despite very large differences in the sample sizes—see S1 Appendix (p. S17). Confidence intervals are in S1 Appendix (p. S12).

https://doi.org/10.1371/journal.pgen.1008612.t002

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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Proportion of narrow-sense chip heritability, A(N) (Eq 37), captured by genome-wide significant SNPs as a function of sample size, N, for phenotypes shown in Figs 1 and Fig 2.

Values for current sample sizes are shown in parentheses. Left-to-right curve order is determined by decreasing . The prediction for education at sample size N = 1.1 million is A(N) = 0.27, so that the proportion of phenotypic variance explained is predicted to be 3.5%, in good agreement with 3.2% reported in [70]. (The curve for AD excludes the APOE locus. For HDL, see S1 Appendix (p. S9) for additional notes).

https://doi.org/10.1371/journal.pgen.1008612.g003

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, ≃ 0.003. In contrast, for AD restricted to chromosome 19, there were an estimated 14 causal SNPs with discoverability ≃ 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 , 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 . For education we found , 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 .

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 doubles, leaving the estimatde number of causal SNPs and heritability as before. For example, for schizophrenia we get π₁ = 5.3 × 10⁻³ and for the reduced reference panel, versus π₁ = 2.8 × 10⁻³ and for the full panel, with heritability remaining essentially the same (37% on the observed scale).

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.