Estimating Genetic Heritability via LDAK Model

We first construct the n × m genotype matrix X, by centering and scaling the allele counts for each SNP according to X_i,j = (S_i,j−2f_j) × [2f_j (1−f_j)]^α/2, where f_j = Σ_iS_i,j/2_n. If w_j and r_j denote the LD weight9 (link) and information score for SNP j, then the LDAK Model for estimating SNP heritability $h_{\text{SNP}}^2={\sigma }_g^2/({\sigma }_g^2+{\sigma }_e^2)$ is: $$\begin{array}{l}{Y}_i={\displaystyle \sum _{k=1}^p{\theta }_k{Z}_{i,k}+}{\displaystyle \sum _{j=1}^m{\beta }_j{X}_{i,j}+e_i,\text{with}}\\ {\beta }_j\sim \text{ℕ}(0,r_j{w}_j{\sigma }_g^2/W),e_i\sim \mathbb{N}(0,{\sigma }_e^2)\\ \text{and}W={\displaystyle \sum _{j=1}^m{r}_j{w}_j{[2f_j(1-f_j)]}^{1+\alpha }.}\end{array}$$
θ_k denotes the fixed-effect coefficient for the kth covariate, β_j and e_i are random-effects indicating the effect size of SNP j and the noise component for Individual i, while ${\sigma }_g^2$ and ${\sigma }_e^2$ are interpreted as genetic and environmental variances, respectively. Note that the introduction of r_j is an addition to the model we proposed in 2012.9 (link)
Model (2) is equivalent to assuming:44 , 45 (link)
$$\mathit{Y}\sim \text{ℕ}(\mathit{Z}\theta ,\mathit{K}{\sigma }_g^2+\mathit{I}{\sigma }_e^2),\text{with}\mathit{K}=\frac{\mathit{X}\mathbf{\Omega }{\mathit{X}}^T}{W},$$ where I is an n × n identity matrix and Ω denotes a diagonal matrix with diagonal entries (r₁w₁, …, r_mw_m). The kinship matrix K, also referred to as a genetic relationship matrix (GRM)1 (link) or genomic similarity matrix (GSM),46 (link) consists of average allelic correlations across the SNPs (adjusted for LD and genotype certainty). Model (3) is typically solved using REstricted Maximum Likelihood (REML), which returns estimates of θ₁, …, θ_p, ${\sigma }_g^2$ and ${\sigma }_e^2.$ 12
The heritability of SNP j can be estimated by $h_j^2={\beta }_j^2\text{Var}(X_j)/\text{Var}(Y),$ which under Model (2), and assuming Hardy-Weinberg Equilibrium,47 (link), 48 has expectation $$\text{𝔼}[h_j^2]=\frac{\text{𝔼}[{\beta }_j^2]\times \text{Var}(X_j)}{\text{Var}(Y)}=\frac{r_j{w}_j{\sigma }_g^2/W\times {[2f_i(1-f_j)]}^{1+\alpha }}{\text{Var}(Y)}.$$ If P₁ and P₂ index two sets of SNPs of size |P₁| and |P₂|, then under the LDAK Model, they are expected to contribute heritability in the ratio W₁ : W₂, where W_l = Σ_j∈Plr_jw_j [2f_j (1−f_j)]^1+α. The GCTA Model corresponds to setting w_j = r_j = 1, in which case W_l = Σ_j∈Pl [2f_j (1−f_j)]^1+α. Most applications of GCTA have further assumed α = −1, so that W_l = |P_l|, which corresponds to the assumption that SNP sets are expected to contribute heritability proportional to the number of SNPs they contain.
Model (2) assumes that all effect-sizes can be described by a single prior distribution. This assumption is relaxed by SNP partitioning. Suppose that the SNPs are divided into tranches P₁, …, P_L of sizes |P₁|, …, |P_L|; typically these will partition the genome, so that each SNP appears in exactly one tranche and Σ_l |P_l| = m, but this is not required. This correspond to generalizing Model (2), so that SNPs in Tranche l have effect-size prior distribution ${\beta }_j\sim \mathbb{N}(0,r_j{w}_j{\sigma }_l^2/W_l).$ Letting $\Sigma ={\sigma }_1^2+\dots +{\sigma }_L^2,$ then $h_{\text{SNP}}^2=\Sigma /(\Sigma +{\sigma }_e^2),$ while ${\sigma }_l^2/\Sigma$ represents the contribution to $h_{\text{SNP}}^2$ of SNPs in Tranche l. This model can equivalently be expressed as $\mathit{\text{Y}}\sim \mathbb{N}(\mathit{\text{Z}}\theta ,{\mathit{\text{K}}}_1{\sigma }_1^2+\dots +{\mathit{\text{K}}}_L{\sigma }_L^2+\mathit{\text{I}}{\sigma }_e^2),$ where K_l represents allele correlations across the SNPs in Tranche l.
For analyses under the LDAK Model, we used LDAK v.5; for analyses under the GCTA Model, we used GCTA v.1.26. For about a third of GCTA-LDMS analyses, the GCTA REML solver failed with the error “information matrix is not invertible,” in which case we rerun using LDAK (while the GCTA and LDAK solvers are both based on Average Information REML,28 (link), 49 (link) subtle differences mean that when using a large number of tranches, one might complete while the other fails). For the few occasions when both solvers failed, we instead used “GCTA-LD” (i.e., SNPs divided only by LD, rather than by LD and MAF), which we found gave very similar results to GCTA-LDMS for traits where both completed (Supplementary Fig. 7). For diseases, we converted estimates of $h_{\text{SNP}}^2$ to the liability scale based on the observed case-control ratio and assumed prevalence.26 (link), 27 (link) In general, we copied the prevalences used by previous studies; however for tuberculosis, where no previous estimate of $h_{\text{SNP}}^2$ is available, we derived an estimate of prevalence from World Health Organization data50 (see Supplementary Note).