Accurate and efficient estimation of local heritability using summary statistics and the linkage disequilibrium matrix

Abstract

Existing SNP-heritability estimators that leverage summary statistics from genome-wide association studies (GWAS) are much less efficient (i.e., have larger standard errors) than the restricted maximum likelihood (REML) estimators which require access to individual-level data. We introduce a new method for local heritability estimation-Heritability Estimation with high Efficiency using LD and association Summary Statistics (HEELS)-that significantly improves the statistical efficiency of summary-statistics-based heritability estimator and attains comparable statistical efficiency as REML (with a relative statistical efficiency >92%). Moreover, we propose representing the empirical LD matrix as the sum of a low-rank matrix and a banded matrix. We show that this way of modeling the LD can not only reduce the storage and memory cost, but also improve the computational efficiency of heritability estimation. We demonstrate the statistical efficiency of HEELS and the advantages of our proposed LD approximation strategies both in simulations and through empirical analyses of the UK Biobank data.

Fig. 1. Comparison of the performance of HEELS with different methods using simulation studies.

Simulated phenotypes using real genotypic data from the UK Biobank, array SNPs on chromosome 22 with MAF > 0.01 (n = 30, 000, p = 9, 205). a Local $h_{SNP}^2$ estimates from HEELS using summary statistics vs GREML using individual-level data. Red-dotted line: y = x. b Analytical SE estimates from HEELS vs GREML. Red-dotted line: y = x. c Distribution of the $h_{SNP}^2$ estimates from 100 simulations using different methods: GREML and BOLT-REML use individual level data; HEELS, GRE and LDSC use summary statistics. Red-dotted line: true $h_{SNP}^2$ of 0.25. The lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles). The upper (lower) whisker extends from the hinge to the largest (smallest) value no further than 1.5 × IQR from the hinge (where IQR is the inter-quartile range). Data beyond the end of the whiskers are called “outlying" points and are plotted individually. d MSEs of the $h_{SNP}^2$ estimates using different methods. e Relative efficiency of different methods compared to GREML.

Fig. 2. Comparison of the performance of different LD approximation strategies.

Simulated phenotypes using real genotypic data from UK Biobank, array SNPs on chromosome 22 with MAF > 0.01 (n = 332, 430, p = 9, 220). Labels of the approximation strategies are explained in Table 2. b denotes the bandwidth of the banded component. r denotes the number of factors in the low-rank component. a Approximation accuracy, measured by $\mid \mid \tilde{\mathbf{R}}\mid {\mid }_F/\mid \mid \mathbf{R}\mid {\mid }_F$, where $\tilde{\mathbf{R}}$ is the approximated LD matrix. Dotted lines are the reference levels: red—85%; blue—95%. b $h_{SNP}^2$ estimates from 100 simulations. Red reference line represents the true heritability value of 0.25. The upper (lower) whisker extends from the mean to the values 1.96 × SE above (below) the mean. The “LR_only" approximations lead to large bias in h² estimation, so are omitted from the comparison.

Fig. 3. Analysis results from application of HEELS to UKB.

a Comparison of heritability estimates and standard errors between BOLT-REML and summary-statistics-based methods (HEELS, GRE, HESS) (n = 332, 340). Each dot represents a local estimate at one locus for a given trait. Red dashed line: y = x. b Heritability enrichment of multiple diseases and traits at putative pleiotropic risk loci. Loci locations on the Y-axis are denoted by the start and end positions (in Mbp). The shade of each box represents the enrichment of heritability for the given region and trait, with log transformation.