Summary statistics

For our real analyses, we use summary statistics from two sets of GWAS. The “25 raw GWAS” (18 binary traits, 7 quantitative, average sample size 9 700; see Supplementary Table 2) are those for which we have access to individual-level data, from either the Wellcome Trust Case Control Consortium¹³ (WTCCC) or the eMERGE Network,¹⁴ and therefore we perform the association analysis ourselves. The “24 summary GWAS” (9 binary traits, 15 quantitative, average sample size 121,000; see Table 1), are those for which we use summary statistics from previously-published analyses.¹⁵

For each of the 25 raw GWAS, we perform strict quality control (see Online Methods), then use REML to estimate how much of the total phenotypic variance explained by SNPs is inflation due to population structure or familial relatedness.^8,16 We estimate that on average 3.1% (range 0.2% to 7.2%) of the variance explained is inflation (Supplementary Fig. 1), indicating that confounding due to population structure or familial relatedness is is modest. For the 24 summary GWAS, we are largely reliant on the quality control decisions of the original authors, which are generally much less strict than ours (Supplementary Table 3). For 4 of the GWAS, imputation info scores are available, so we exclude SNPs with score < 0.95; for the remaining 20 GWAS, we instead restrict to the 4,555,718 SNPs present in the eMERGE data, as we observe that these SNPs tend to be well-imputed in Caucasian GWAS (Supplementary Table 4).

The importance of the heritability model

We begin by analyzing summary statistics from simulated phenotypes. Using the eMERGE data^14,17 (25,875 individuals recorded for 4,555,718 SNPs), we generate 200 phenotypes each with 2,000 causal SNPs, h²_SNP = 0.5, no confounding bias, and such that consecutive pairs (Phenotypes 1 & 2, 3 & 4, etc) have genetic correlation 0.5. For Phenotypes 1-100, we sample causal SNP effect sizes according to the GCTA Model, for Phenotypes 101-200, according to the LDAK Model. We perform single-SNP analysis for each phenotype, then analyze the resulting summary statistics using LDSC-Zero, LDSC, SumHer-Zero and SumHer-GC. Selected results are provided in Figure 1, with full results in Supplementary Figures 2 & 3.

Open in a separate window

Figure 1

Importance of the heritability model.

We generated 100 phenotypes assuming the GCTA heritability model, 100 assuming the LDAK heritability model, then analyzed each using LDSC-Zero, LDSC, SumHer-Zero and SumHer-GC (see main text for details of the heritability models and methods). a, Average estimates of h²_SNP (true h²_SNP is 0.5). b, Average estimates of the enrichment of heritability in conserved regions (true enrichment is 1). c, Average estimates of genetic correlation between pairs of phenotypes (true correlation is 0.5). In all plots, vertical line segments mark 95% confidence intervals for the average estimates.

First we examine results from analyzing each phenotype individually. Figure 1a shows that, as expected, accurate estimates of h²_SNP are returned when phenotypes are analyzed assuming the matching heritability model (i.e., when GCTA phenotypes are analyzed using LDSC-Zero or LDSC and LDAK phenotypes are analyzed using SumHer-Zero or SumHer-GC), but that using a different heritability model can result in very poor estimates; SumHer-Zero tends to over-estimate h²_SNP by about 20% when applied to GCTA phenotypes, while LDSC-Zero tends to under-estimate h²_SNP by about 40% when applied to LDAK phenotypes. As shown in Supplementary Figure 2a, LDSC infers that there is only slight confounding when used on GCTA phenotypes (average intercept 1.008, s.d. 0.001), and therefore its estimates of h²_SNP closely match those from LDSC-Zero. However, when used on LDAK phenotypes, LDSC wrongly infers that much of the causal signal is in fact confounding (average intercept 1.096, s.d. 0.001), and as a result, its estimates of h²_SNP are on average about 60% lower than those from LDSC-Zero and about 75% lower than the true value.

Figure 1b reports the estimated enrichment of SNPs in conserved regions. As causal SNPs were picked at random from across the genome, the true enrichment is one. Again, we see that assuming the correct heritability model produces reliable estimates, but assuming the wrong model can result in misleading conclusions. In particular, we find that when LDSC-Zero or LDSC are used to analyze LDAK phenotypes, they infer that conserved regions are over 2-fold enriched for heritability. We chose to focus on conserved regions as this was the category that, by applying LDSC to real data, Finucane et al.² found to be most enriched. These simulations show that a substantial portion of the enrichment observed by Finucane et al. could be an artifact of misspecifying the heritability model.

Figure 1c compares results from analyzing consecutive pairs of phenotypes. Whereas LDSC-Zero and SumHer-Zero only produce accurate estimates of genetic correlation when applied to GCTA and LDAK phenotypes, respectively, we find that LDSC and SumHer-GC produce accurate estimates for both sets of phenotypes (Supplementary Fig. 2f shows why this is the case). Even so, highest precision is achieved when the correct heritability model is assumed; the s.d. of SumHer-GC estimates is on average about 50% higher than that of LDSC estimates when analyzing GCTA phenotypes, but about 50% lower when analyzing LDAK phenotypes.

Determining the most appropriate heritability model

With access to individual-level data, we can compare how well different heritability models fit real data using the likelihood from REML analysis.⁵ With only summary statistics, we can instead use the likelihood from SumHer. Supplementary Table 5 shows that across the 25 raw GWAS, the SumHer log likelihood is on average 17 nats higher under the LDAK Model than under the GCTA Model (for comparison, the average difference in REML log likelihood is 19 nats; Supplementary Table 6). Across the 24 summary GWAS, the SumHer log likelihood is on average 76 nats higher under the LDAK Model (Supplementary Table 7).

Rather than comparing based on likelihoods, an easier-to-visualize approach is to run SumHer using a hybrid heritability model containing both the GCTA and LDAK Models, and allow the data to decide the relative weighting of each. Specifically, we use Hybrid-Zero (or Hybrid-GC if the test statistics are confounded), with the focus on estimating p, the “proportion of LDAK”, in Model (2). Figure 2a demonstrates proof of principle using our simulated phenotypes: we see that Hybrid-Zero correctly estimates p ≈ 0 when applied to the GCTA phenotypes, p ≈ 1 for the LDAK phenotypes, and p ≈ 0.5 for “Hybrid Phenotypes” (created by summing Phenotypes 1 & 101, 2 & 102, etc.). For Figure 2b and Supplementary Table 5, we apply Hybrid-Zero to the 25 raw GWAS; the average estimate of p is 1.03 (s.d. 0.02), indicating that the data overwhelmingly support the LDAK Model over the GCTA Model. For Figure 2c and Supplementary Table 7, we apply Hybrid-GC to the 24 summary GWAS; the average estimate of p is 0.85 (s.d. 0.01), again strongly supporting the LDAK Model (in the Discussion, we consider why this estimate is lower than that for the 25 raw GWAS).

Open in a separate window

Figure 2

Comparing the GCTA and LDAK heritability models.

These analyses use Hybrid-Zero and Hybrid-GC, versions of SumHer that assign weights 1-p and p to the GCTA and LDAK heritability models, respectively. a, Average estimates of p from Hybrid-Zero for GCTA phenotypes (true p = 0), LDAK phenotypes (true p = 1) and hybrid phenotypes (true p = 0.5). b, Estimates of p from Hybrid-Zero for the 25 raw GWAS. Colors distinguish between the 13 WTCCC, 5 binary eMERGE and 7 quantitative eMERGE traits (black denotes the 25-trait average). A precise estimate of p was not possible for shingles (Segment 17), due to the trait having very low h²_SNP. c, Estimates of p from Hybrid-GC for the 24 summary GWAS. Colors distinguish between the 9 binary and 15 quantitative traits (black denotes the 24-trait average). In all plots, vertical line segments mark 95% confidence intervals.

Population structure and relatedness

Supplementary Tables 2, 8 & 9 reports estimates of confounding bias, SNP heritability, enrichments and genetic correlations for the 25 raw GWAS. As a consequence of our strict quality control, SumHer-GC finds only slight confounding bias; its average estimate of the scaling factor C is 1.005 (s.d. 0.002). We now construct GWAS with substantial confounding. First, for each of the 13 WTCCC GWAS, we introduce population structure by replacing 2000 of the controls (on average 54%) with 2000 individuals from POBI¹⁸ (People of the British Isles); this generates population structure because, although both WTCCC and POBI individuals were recruited from the UK, the latter predominately came from isolated, rural regions (Supplementary Fig. 4). Second, for each of the 12 eMERGE GWAS, we no longer filter individuals based on relatedness; this leads to the retention of approximately 1650 pairs of relatives (Supplementary Fig. 5). Supplementary Tables 10 & 11 confirm that in both cases, SumHer-GC now detects significant confounding; its average estimates of C are 1.022 (s.d. 0.003) and 1.020 (s.d. 0.003), respectively.

Genomic control and mixed-model association analysis

Although its use is declining, it remains that the majority of published GWAS have performed genomic control (divided test statistics by the GIF) at least once in their analyses.¹¹ As the GIF tends to over-estimate confounding bias,¹⁹ genomic control will tend to result in deflated test statistics. Supplementary Figures 6 & 7 show that, if not accounted for, genomic control will result in under-estimation of h²_SNP and inaccurate estimates of heritability enrichments, but that reliable estimates can be obtained by allowing for multiplicative inflation of test statistics. SumHer, like LDSC, is designed to be used with results from classical (least-squared) regression. However, with the development of software such as Fast-LMM, GCTA-LOCO and Bolt-LMM^20–22 it is now common for GWAS to instead use mixed-model association analysis.^23,24 Supplementary Figures 6, 8 & 9 show that when estimating h²_SNP and enrichments, the impact of mixed-model analysis is similar to, albeit less severe than, that of genomic control. However, as with genomic control, reliable estimates can be obtained by allowing for multiplicative inflation of test statistics.

Estimates of confounding bias and SNP heritability

Table 1 reports estimates of confounding bias for each of the 24 summary GWAS. LDSC estimates the average intercept to be 1.042 (s.d. 0.002), indicating that the test statistics tend to be inflated. By contrast, SumHer-GC estimates the average scaling factor to be 0.929 (s.d. 0.003), indicating that the test statistics tend to be deflated. The widespread deflation observed by SumHer-GC is consistent with the fact that at least 15 of the 24 summary GWAS performed genomic control (Supplementary Table 3); in particular, we note that the four GWAS estimated to have lowest C (0.55 for body mass index, 0.68 for HDL cholesterol, 0.73 for LDL cholesterol and 0.70 for triglyceride levels), are all meta-analyses that used genomic control both before and after combining results across cohorts.^25,26 Table 1 also reports the number of independent loci with P < 5×10-8. Without adjustment, there are on average 86 loci per trait. If we adjust based on the results of LDSC (i.e., for each trait divide test statistics by the LDSC estimate of confounding bias), the average reduces to 62, whereas if we adjust based on the results of SumHer-GC, it increases to 118. For these counts, we defined two significant SNPs as dependent if they are within 1 cM and have r²_jl > 0.1, but we find that results are very similar if instead we increase the window size to 3 cM or use r²_jl > 0.2 (Supplementary Table 12).

Table 1 shows that across the 24 traits, estimates of h²_SNP from SumHer-GC are on average 2.7 (s.d. 0.05) times higher than those from LDSC. As shown in Supplementary Table 13, this difference is primarily due to changing the heritability model, rather than changing the confounding model; for example, were we to switch from the GCTA to LDAK Model but to continue using an additive model for confounding inflation, estimates would still be on average 2.3 (s.d. 0.04) times higher.

Estimates of functional enrichments and genetic correlations

Figure 3a and Supplementary Table 14 report estimates of enrichments for the 24 functional categories, averaged across the 24 traits. We see striking differences between the estimates from LDSC and SumHer-GC. For example, LDSC estimates of enrichments range from -1.1 to 9.4, whereas SumHer-GC estimates range from 0.78 to 2.0. The estimated enrichment of a category is its estimated share of h²_SNP divided by its expected share under the assumed heritability model.²⁷ Supplementary Table 14 shows that changing from the GCTA to LDAK model has a small impact on the expected share of each category (the denominator), but typically has a large impact on the estimated share (the numerator). Supplementary Figure 10 and Supplementary Tables 15 & 16 show that the large differences between LDSC and SumHer-GC estimates of enrichments remain if we calculate the former using the LDSC software¹ (in this setting, LDSC is often referred to as stratified LDSC or S-LDSC) rather than using our implementation of LDSC in SumHer, and likewise if we vary the SNP sets (the authors of LDSC recommend that the reference panel contains as many SNPs as possible, but that only HapMap3 SNPs²⁸ with MAF ≥ 0.05 are used when performing the regression^1,2), or if we use the 75-part model proposed by Gazal et al.²⁹

Open in a separate window

Figure 3

Functional enrichments across the 24 summary GWAS.

a, Average estimates of enrichments for the 24 functional categories from LDSC (red bars) and from SumHer-GC (blue bars). b, Average estimates of enrichments from SumHer-GC, based either on the 9 binary or on the 15 quantitative traits. In both plots, horizontal and vertical line segments mark 95% confidence intervals.

Based on the results from SumHer-GC, we conclude that conserved regions^30,31 and transcription start sites³² are most enriched for heritability; they are estimated to contribute 7.1% (s.d. 0.2) and 3.6% (s.d. 0.2) of h²_SNP, respectively, 1.95 (s.d. 0.07) and 1.97 (s.d. 0.09) times higher than their expected contributions under the LDAK Model. Repressed region are the only category significantly depleted; although they are estimated to contribute 35% (s.d. 0.5) of h²_SNP, this is only 0.78 (s.d. 0.01) times their expected contribution.

Supplementary Figure 11 and Supplementary Table 17 provide estimates of genetic correlation for the 276 pairs of traits. As predicted by our simulation study, there is strong concordance between estimates from LDSC and SumHer-GC, but the SumHer-GC estimates are more precise; for example, across the 38 pairs of traits that both methods find to be significantly correlated (P < 0.05/276), the s.d. of SumHer-GC estimates is on average about 20% lower than that of LDSC estimates.

The heritability model matters

The heritability model describes the prior assumptions regarding how heritability is distributed across the genome. We began by showing that estimates of confounding bias, h²_SNP and enrichments are sensitive to the assumed heritability model, and that using an unrealistic heritability model can result in misleading conclusions. Next we evaluated the GCTA and LDAK heritability models on real data. Not only did we demonstrate that the LDAK Model is substantially more realistic than the GCTA Model, but also our prediction analysis showed that the LDAK Model is sufficiently realistic to be useful. Nonetheless, it is important to recognize that the LDAK heritability model is not perfect, and as such SumHer provides tools for testing and comparing different heritability models on large-scale GWAS data.

Hybrid heritability models

As well as comparing the GCTA and LDAK Models based on likelihood, we also ran SumHer using a hybrid heritability model where the fractions 1-p and p indicate the proportions of GCTA and LDAK, respectively; across the 25 raw GWAS, the average estimate of p was 1.03 (s.d. 0.02), while across the 24 summary GWAS, the average estimate of p was 0.85 (s.d. 0.01). Although both estimates of p are greater than 0.5, consistent with the LDAK Model being more realistic than the GCTA Model, the latter suggests the hybrid model might improve on the LDAK Model. We recommend using only the LDAK Model for two reasons. Firstly, estimates from the LDAK Model tend to be more precise than those from the hybrid model. This is particularly true when estimating enrichments; for example, Supplementary Table 19 shows that across the 24 summary GWAS, the s.d. of SumHer-GC estimates is on average about 70% lower than that of Hybrid-GC estimates. Secondly, Supplementary Figure 12 shows that introducing poorly-genotyped SNPs makes traits appear more GCTA-like, indicating that the hybrid model is more sensitive to genotyping errors. In any case, Supplementary Figure 13 shows that even if we had instead preferred results from the hybrid model, it would not have affected our overall conclusions.

Confounding bias in GWAS

In addition to examining the choice of heritability model, we also considered how to model inflation of test statistics due to confounding. Whereas LDSC uses an additive model for confounding inflation, we prefer a multiplicative model, due to the fact that genomic control (which scales test statistics multiplicatively) has been widely used by published GWAS (we provide further support for our choice in Supplementary Fig. 14 and Supplementary Table 13). We recognize that as genomic control becomes less common, it will be necessary to check whether the multiplicative model remains preferable to the additive model.

Comparing heritability models

The (weighted) natural log likelihood is

To evaluate L(S | h²_SNP, D) requires values for h²_SNP and D. While the natural choice is to use the values after the final iteration, in order to compare different heritability models based on L(S | h²_SNP, D), we must use the same D for each (else models leading to lower D_jj would have an unfair advantage). Therefore, SumHer reports $L(S|\widehat{h_{SNP}^2,D^0)},$ where $\widehat{h_{SNP}^2}$ is the final estimate of h²_SNP, but D⁰ is the initial weight matrix (obtained by setting 1 / D_jj = ∑_lNj r²_jl); Supplementary Tables 5 & 6 show that for the 25 raw GWAS, comparisons of heritability models based on $L(S|\widehat{h_{SNP}^2,D^0)}$ align closely with those based on REML likelihood.

Estimating confounding bias

When assuming that confounding inflates test statistics multiplicatively, we replace Equation (5) with E[S_j] ≈ C(1 + u_j h²_SNP), then estimate C and h²_SNP by jointly regressing S_j on 1 and u_j. To instead copy LDSC and assume confounding causes additive inflation, we replace Equation (5) with E[Sj] ≈ 1 + n_j a + u_j h²_SNP, then estimate a and h²_SNP by jointly regressing (S_j - 1) on n_j and u_j (note that in the main text we assumed n_j = n, allowing us to replace n_j a by A).

Estimating enrichments

Suppose we have K categories; let I_jk {0,1} indicate whether SNP j belongs to Category k. We wish to estimate h²_{Cat k} = ∑ _j I_jkh²_j, the heritability contributed by SNPs in Category k. We now use the heritability model

where q_j I_jk h²_{Part k} / Q_k represents the contribution of Category k to the expected heritability of SNP j, and h²_SNP = h²_{Part 1} + … + h²_{Part K}. We consider two scenarios: in Scenario 1, the K categories partition the genome (∑_k I_jk = 1); in Scenario 2, the first K-1 categories correspond to annotations, and the Kth is the base category containing all SNPs (I_jK = 1, Q_K = m). Using Equation (6) ensures that when there is no enrichment, E[h²_j] = q_j h²_SNP / Q for both scenarios. To appreciate why, consider that for Scenario 1, h²_{Part k} = h²_{Cat k}, so that when no categories are enriched, h²_{Part k} = Q_k h²_SNP / Q; for Scenario 2, h²_{Part 1}, … , h²_{Part K-1} indicate how much each annotation increases the SNP heritability from h²_{Part K}, so that when there is no enrichment, h²_{Part 1} = … = h²_{Part K-1} = 0 and h²_{Part K} = h²_SNP. Now when we replace h²_j in Equation (4) by its expected value, we obtain

and therefore we can estimate h²_{Part 1}, …, h²_{Part K} by jointly regressing (S_j - 1) on u_j1, …, u_jK. Given these, our estimate of h²_{Cat k} is ∑_j,k′ (I_jk q_j I_jk′ h²_{Part k’} / Q_k’), which we then divide by Q_k h²_SNP / Q to get an estimate of the enrichment of Category k.

Estimating genetic correlation

Suppose we have summary statistics from two GWAS. Instead of χ²(1) test statistics, we now use (signed) Z-statistics. Let Z_Aj and Z_Bj denote the two Z-statistics for SNP j, computed using n_Aj and n_Bj individuals, respectively, of which n_Cj were common to both GWAS (if the two GWAS were independent, n_Cj = 0). We assume

where c_AB is the phenotypic correlation between the two traits and h²_AB is their genetic covariance. This equation matches that used by LDSC,³ except that we have replaced r²_jl / m by q_j r²_jl / Q. By regressing Z_AjZ_Bj on u′_j, we obtain an estimate of h²_AB, which we then divide by estimates of $\sqrt{h_{SNP}^2}$ for each trait to get an estimate of their genetic correlation.

Impact of changing the confounding model

Note that when per-SNP sample sizes are constant, n_j = n, so we can write 1 + n_j a + u_j h²_SNP as C (1+ u_j h²_SNP / C), where C = 1 + n a = 1 + A. Therefore, given the heritability model (which determines the u_j), the estimate of confounding bias will be the same whether we assume additive inflation of test statistics (and estimate 1+A) or assume multiplicative inflation (and estimate C), but estimates of h²_SNP from the additive model will be C times those from the multiplicative model (estimates of enrichments and genetic correlations will be unaffected, as these are ratios of heritabilities so the factor C cancels).

Reference panel

Our primary reference panel is 404 non-Finnish, European individuals from the 1000 Genomes Project.¹² In Supplementary Figure 16, we repeat our analyses of the 24 summary GWAS using instead 8,850 unrelated Caucasian individuals from the Health & Retirement Study. The results are almost identical, indicating that despite its small size, it suffices to construct reference panels from the 1000 Genomes Project data.

Regions of extreme LD and large-effect loci

Estimates can be unduly affected by regions of extreme LD and by SNPs with disproportionately large effect size.^1,4,5 Therefore, for all analyses, we exclude SNPs within the MHC (Chromosome 6: 25-34 Mb), as well as SNPs which individually explain >1% of phenotypic variation (S_j > n_j / 99), and SNPs in LD with these (within 1 cM and r²_jl > 0.1); the latter resulted in the exclusion of SNPs for 10 of the 25 raw GWAS and for 4 of the 24 summary GWAS (Supplementary Table 20). Excluding large-effect loci will result in under-estimation of h²_SNP; this has little consequence for our analyses, as we are primarily interested in comparing models, but if preferred can be avoided by re-introducing the large-effect SNPs as fixed effects.⁵

Binary phenotypes

So far, we have implicitly assumed the trait is quantitative. If instead it is binary (i.e., for case-control GWAS), there are two considerations. Firstly, heritability estimates now correspond to the observed scale; SumHer will convert them to the liability scale if provided with the prevalence and case-control ratio.^51,52 Secondly, the p-values may have come from logistic regression instead of linear regression; Supplementary Figure 17 shows that because linear and logistic p-values closely match for SNPs with small or moderate effect, this has limited impact.

Polygenic risk scoring

To predict phenotypes, we construct PRS of the form ∑_j β_j X′_j, where the vector X′_j contains standardized genotypes for SNP j (obtained by centering X_j, then scaling to have variance 1). Without loss of generality, we assume phenotypes have also been standardized, in which case the estimate of β_j from classical linear regression is ρ_j, the correlation observed between SNP j and the phenotype, and has variance (1 – ρ²_j) / n_j. We can infer ρ_j from the Wald test Z-statistic $p_j/\sqrt{(1-p_j^2)/n_j}.$ For the Classical PRS, our estimate of β_j is ρ_j. For the Bayesian PRS, we use the prior distribution β_j ~ (0, σ²_j), where σ²_j = (q_j + ∑_lNj q_l r²_jl) h²_SNP / Q; this choice is motivated by recognizing that for the “true PRS”, β²_j = v²_j, the heritability tagged by SNP j. We then estimate β_j by its posterior mean, a shrunken version of the classical estimate. To calculate this, we approximate the likelihood by (ρ_j, (1 – ρ²_j) / n_j), then the posterior mean equals ρ_j σ²_j / (σ²_j + (1 – ρ²_j) / n_j).^53,54 There are two technical issues. Firstly, to construct each Bayesian PRS requires a value for h²_SNP, so we used the corresponding estimates from LDSC-Zero and SumHer-Zero; Supplementary Table 18 confirms that predictions are approximately the same if we instead use the estimates from LDSC and SumHer-GC, or agnostically set h²_SNP = 0.5. Secondly, for the Bayesian PRS, we performed clumping (identified pairs of SNPs within 1 cM with r²_jl > 0.5, then discarded the least significant), whereas for the Classical PRS we did not; this was because we found that clumping benefited all Bayesian PRS, but did not benefit the Classical PRS.

Choosing a heritability model

Although the heritability model is defined in terms of h²_j, as is clear from Equation (4), its fit depends only on how well it models v²_j = h²_j + ∑_lNj r²_jl h²_l. Therefore, when specifying q_j, the focus should be on accurately describing how v²_j varies across the genome. LDSC assumes the GCTA Model,^5,7 (q_j = 1), which implies that v²_j correlates with local levels of LD. Despite its widespread use in statistical genetics (e.g., it is implicitly assumed by any regression method where SNPs are standardized and the same penalty function / prior distribution is applied to each), the GCTA Model poorly reflects real data.⁵ Even for traits where a correlation is observed between v²_j and S_j (those with p < 1 in Figures 2b & 2c), the correlation is substantially weaker than predicted by the GCTA Model (hence p >> 0), and partly a consequence of including lower-certainty SNPs in the GWAS (Supplementary Fig. 12). We recommend using the LDAK Model. Originally, this took the form q_j = w_j, where the weightings w_j are calculated based on the assumption that E[v²_j] is constant.⁶ We recently updated it to q_j = w_j [f_j(1 - f_j)]^0.75 r_j, reflecting that v²_j tends to depend on MAF, f_j, and genotype certainty, r_j (we did not use the latter as either r_j ≈ 1 or was unknown).⁵ The exponent 0.75 was chosen as this value performed well in a previous analysis of 42 human traits. In Supplementary Figure 18, we confirm 0.75 also performs well for both the 25 raw and 24 summary GWAS.

Quality control

To prepare the 13 WTCCC GWAS, we used our previously-described protocol.⁵ In summary, after excluding apparent population outliers, individuals with extreme missingness or heterozygosity and SNPs with MAF < 0.01, call rate < 0.95 or Hardy-Weinberg P < 1 × 10^-6, we phased using SHAPEIT, then imputed using IMPUTE2 with the 1000 Genomes Project Phase 3 (2014) reference panel.^12,55,56 When merging case and control datasets, we converted genotype probabilities to hard calls using a threshold of 0.95, then retained only autosomal SNPs that in all cohorts had MAF ≥ 0.01 and info score ≥ 0.99. Finally, we thinned individuals, so no pair remained with allelic correlation⁵⁷ > 0.05. For the 12 eMERGE GWAS, data were provided post-imputation; this was also performed using SHAPEIT and IMPUTE2, but used the 1000 Genomes Project Phase 2 reference panel.¹⁷ We converted genotype probabilities to hard calls using a certainty threshold of 0.95, then retained only biallelic SNPs with MAF ≥ 0.01, call rate ≥ 0.95, info score ≥ 0.95 and whose genomic positions matched those in the 1000 Genomes Project. Finally, we excluded individuals whose genotypes were inconsistent with those of non-Finnish Europeans from the 1000 Genomes Project (Supplementary Fig. 19), and those whose ethnicity was reported as “Hispanic or Latino”, then filtered until no pair remained with allelic correlation⁵⁷ > 0.05 (which left 25,875 individuals). For the 24 summary GWAS, quality control varied by study (see Table 1 for references), but was in general much less strict than we would have performed had we access to the individual-level data¹⁶ (for example, the most common info score thresholds were 0.3 and 0.5). Info scores were only available for Crohn’s Disease, inflammatory bowel disease, schizophrenia and ulcerative colitis; for these traits we restricted to SNPs with score ≥ 0.95, for the remainder, we restricted to SNPs with info score ≥ 0.95 in the eMERGE data (Supplementary Table 4).

Association analysis

For the 25 raw GWAS, we performed the association analysis using linear regression (regardless of whether the trait was quantitative or binary), including as covariates sex and ten principal components (five derived from the reference panel, five from the 1000 Genomes Project). The 24 summary GWAS varied in whether they used classical or mixed-model regression, and whether they performed a meta or mega analysis (Supplementary Table 3).

Software

SumHer is implemented in Version 5.0 of the LDAK software. In Supplementary Figure 20 we confirm that LDSC estimates computed using SumHer match those computed using Version 1.0.0 of the LDSC Software.¹

We thank A Price, H Finucane, P O’Reilly and M Speed for helpful discussions. Access to Wellcome Trust Case Control Consortium data was authorized as work related to the project “Genome-wide association study of susceptibility and clinical phenotypes in epilepsy”, access to eMERGE Network data was granted under dbGaP Project 14422, “Comprehensive testing of SNP-based prediction models”, while access to the Health and Retirement Study was granted under dbGaP Project 15139, “Developing summary-statistic tools for analysing genetic association study data”. D.S. is funded by the UK Medical Research Council under grant MR/L012561/1, by the European Unions Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement number 754513, by Aarhus University Research Foundation (AUFF) and the Independent Research Fund Denmark under Project 7025-00094B. The eMERGE Network was initiated and funded by NHGRI through the following grants: U01HG006828 (Cincinnati Childrens Hospital Medical Center/Boston Childrens Hospital); U01HG006830 (Childrens Hospital of Philadelphia); U01HG006389 (Essentia Institute of Rural Health, Marshfield Clinic Research Foundation and Pennsylvania State University); U01HG006382 (Geisinger Clinic); U01HG006375 (Group Health Cooperative); U01HG006379 (Mayo Clinic); U01HG006380 (Icahn School of Medicine at Mount Sinai); U01HG006388 (Northwestern University); U01HG006378 (Vanderbilt University Medical Center); and U01HG006385 (Vanderbilt University Medical Center serving as the Coordinating Center). The Health and Retirement Study genetic data is sponsored by the National Institute on Aging (grant numbers U01AG009740, RC2AG036495, and RC4AG039029) and was conducted by the University of Michigan. Analyses were performed with the use of the UCL Computer Science Cluster and the help of the CS Technical Support Group, as well as the use of the UCL Legion High-Performance Computing Facility (Legion@UCL) and associated support services.