Improved heritability estimation from genome-wide SNPs

Abstract

Estimation of narrow-sense heritability, h(2), from genome-wide SNPs genotyped in unrelated individuals has recently attracted interest and offers several advantages over traditional pedigree-based methods. With the use of this approach, it has been estimated that over half the heritability of human height can be attributed to the ~300,000 SNPs on a genome-wide genotyping array. In comparison, only 5%-10% can be explained by SNPs reaching genome-wide significance. We investigated via simulation the validity of several key assumptions underpinning the mixed-model analysis used in SNP-based h(2) estimation. Although we found that the method is reasonably robust to violations of four key assumptions, it can be highly sensitive to uneven linkage disequilibrium (LD) between SNPs: contributions to h(2) are overestimated from causal variants in regions of high LD and are underestimated in regions of low LD. The overall direction of the bias can be up or down depending on the genetic architecture of the trait, but it can be substantial in realistic scenarios. We propose a modified kinship matrix in which SNPs are weighted according to local LD. We show that this correction greatly reduces the bias and increases the precision of h(2) estimates. We demonstrate the impact of our method on the first seven diseases studied by the Wellcome Trust Case Control Consortium. Our LD adjustment revises downward the h(2) estimate for immune-related diseases, as expected because of high LD in the major-histocompatibility region, but increases it for some nonimmune diseases. To calculate our revised kinship matrix, we developed LDAK, software for computing LD-adjusted kinships.

Figure 1

Investigation of the Robustness of ${\hat{h}}^2$ to Assumptions of Polygeneity (A) The distribution of ${\hat{h}}^2$ for different numbers of causal variants, from one up to “ALL” (all 81,327 SNPs), with the use of the standard kinship matrix A (left) and the weighted kinship matrix A^∗ (right). Boxes indicate interquartile ranges, colors correspond to simulated h² (red, 0.5; green, 0.8), and whiskers span the full range except for outliers, indicated with circles. (B) The layout matches that of (A), but now the boxes correspond to the REML SD estimates calculated by GCTA, and the purple lines mark the empirical SD estimates based on the 50 replicates.

Figure 2

Investigation of the Robustness of ${\hat{h}}^2$ to Assumptions of the Relationship between Effect-Size Variance and MAF Phenotypes were simulated with each of four models (indexed by α₁) for the relationship between effect-size variance and MAF (Equation 5). Analysis was performed with each of the same four models (indexed by α₂) when allele counts were standardized. Boxes indicate interquartile ranges of ${\hat{h}}^2$. Colors correspond to simulated h² (red, 0.5; green, 0.8), and gray boxes indicate that the analysis model matches the simulation model (α₁ = α₂).

Figure 3

Distributions of ${\hat{h}}^2$ with and without Adjustment for LD The x axis indicates the relative levels of tagging of the causal variants. The boxes indicate interquartile ranges of ${\hat{h}}^2$ under SNP-based mixed-model analysis using A (left) or A^∗ (right). Colors correspond to simulated h² (red, 0.5; green, 0.8), and gray boxes indicate that causal variants were chosen at random without regard to tagging.