Discovering nonadditive heritability using additive GWAS summary statistics
Abstract
LD score regression (LDSC) is a method to estimate narrowsense heritability from genomewide association study (GWAS) summary statistics alone, making it a fast and popular approach. In this work, we present interactionLD score (iLDSC) regression: an extension of the original LDSC framework that accounts for interactions between genetic variants. By studying a wide range of generative models in simulations, and by reanalyzing 25 wellstudied quantitative phenotypes from 349,468 individuals in the UK Biobank and up to 159,095 individuals in BioBank Japan, we show that the inclusion of a cisinteraction score (i.e. interactions between a focal variant and proximal variants) recovers genetic variance that is not captured by LDSC. For each of the 25 traits analyzed in the UK Biobank and BioBank Japan, iLDSC detects additional variation contributed by genetic interactions. The iLDSC software and its application to these biobanks represent a step towards resolving further genetic contributions of sources of nonadditive genetic effects to complex trait variation.
Editor's evaluation
This study provides a valuable investigation into whether phenotypic variance due to interactions between genetic variants can be measured using genomewide association summary statistics. The authors present a convincing method, iLDSC, that uses statistics on the correlations between genotypes at different loci (linkage disequilibrium) to estimate the phenotypic variance explained by both additive genetic effects and pairwise interactions.
https://doi.org/10.7554/eLife.90459.sa0Introduction
Heritability is defined as the proportion of phenotypic trait variation that can be explained by genetic effects (BulikSullivan et al., 2015b, BulikSullivan et al., 2015a, Shi et al., 2016). Until recently, studies of heritability in humans have been reliant on typically small sized family studies with known relatedness structures among individuals (Zaitlen et al., 2013; Polderman et al., 2015). Due to advances in genomic sequencing and the steady development of statistical tools, it is now possible to obtain reliable heritability estimates from biobankscale data sets of unrelated individuals (BulikSullivan et al., 2015b; Shi et al., 2016; Hou et al., 2019; Pazokitoroudi et al., 2020). Computational and privacy considerations with genomewide association studies (GWAS) in these larger cohorts have motivated a recent trend to estimate heritability using summary statistics (i.e. estimated effect sizes and their corresponding standard errors). In the GWAS framework, additive effect sizes and standard errors for individual single nucleotide polymorphisms (SNPs) are estimated by regressing phenotype measurements onto the allele counts of each SNP independently. Through the application of this approach over the last two decades, it has become clear that many traits have a complex and polygenic basis—that is, hundreds to thousands of individual genetic loci across the genome often contribute to the genetic basis of variation in a single trait (Yengo et al., 2018).
Many statistical methods have been developed to improve the estimation of heritability from GWAS summary statistics (BulikSullivan et al., 2015b, Shi et al., 2016, Speed and Balding, 2019, Song et al., 2022). The most widely used of these approaches is linkage disequilibrium (LD) score regression and the corresponding LDSC software (BulikSullivan et al., 2015b), which corrects for inflation in GWAS summary statistics by modeling the relationship between the variance of SNPlevel effect sizes and the sum of correlation coefficients between focal SNPs and their genomic neighbors (i.e. the LD score of each variant). The formulation of the LDSC framework relies on the fact that the expected relationship between chisquare test statistics (i.e. the squared magnitude of GWAS allelic effect estimates) and LD scores holds when complex traits are generated under the infinitesimal (or polygenic) model which assumes: (i) all causal variants have the same expected contribution to phenotypic variation and (ii) causal variants are uniformly distributed along the genome. Initial simulations in BulikSullivan et al. showed that violations of these assumptions can be tolerated to a point, but begin to affect the estimation of narrowsense heritability once a certain proportion of variants have nonzero effects. Importantly, the estimand of the LDSC model is the proportion of phenotypic variance attributable to additive effects of genotyped SNPs. The main motivation behind the LDSC model is that, for polygenic traits, many marker SNPs tag nonzero effects. This may simply arise because some of these SNPS are in LD with causal variants (BulikSullivan et al., 2015b) or because their statistical association is the product of a confounding factor such as population stratification.
As of late, there have been many efforts to build upon and improve the LDSC framework. For example, recent work has shown that it is possible to estimate the proportion of phenotypic variation explained by dominance effects (Palmer et al., 2023) and local ancestry (Chan et al., 2023) using extensions of the LDSC model. One limitation of LDSC is that, in practice, it only uses the diagonal elements of the squared LD matrix in its formulation which, while computationally efficient, does not account for information about trait architecture that is captured by the offdiagonal elements. This tradeoff helps LDSC to scale genomewide, but it has also been shown to lead to heritability estimates with large standard error (Ning et al., 2020, Zhang et al., 2021, Song et al., 2022). Recently, newer approaches have attempted to reformulate the LDSC model by using the eigenvalues of the LD matrix to leverage more of the information present in the correlation structure between SNPs (Shi et al., 2016, Song et al., 2022).
In this paper, we show that the LDSC framework can be extended to estimate greater proportions of genetic variance in complex traits (i.e. beyond the variance that is attributable to additive effects) when a subset of causal variants is involved in a genebygene (G×G) interaction. Indeed, recent association mapping studies have shown that G×G interactions can drive heterogeneity of causal variant effect sizes (Patel et al., 2022). Importantly, nonadditive genetic effects have been proposed as one of the main factors that explains ‘missing’ heritability—the proportion of heritability not explained by the additive effects of variants (Eichler et al., 2010).
The key insight we highlight in this manuscript is that SNPlevel GWAS summary statistics can provide evidence of nonadditive genetic effects contributing to trait architecture if there is a nonzero correlation between individuallevel genotypes and their statistical interactions. We present the ‘interactionLD score’ regression model or iLDSC: an extension of the LDSC framework which recovers ‘missing’ heritability by leveraging this ‘tagged’ relationship between linear and nonlinear genetic effects. To validate the performance of iLDSC in simulation studies, we focus on synthetic trait architectures that have been generated with contributions stemming from secondorder and cisacting statistical SNPbySNP interaction effects; however, note that the general concept underlying iLDSC can easily be extended to other sources of nonadditive genetic effects (e.g. genebyenvironment interactions). The main difference between iLDSC and LDSC is that the iLDSC model includes an additional set of ‘cisinteraction’ LD scores in its regression model. These scores measure the amount of phenoytpic variation contributed by genetic interactions that can be explained by additive effects. In practice, these additional scores are efficient to compute and require nothing more than access to a representative pairwise LD map, same as the input required for LD score regression.
Through extensive simulations, we show that iLDSC recovers substantial nonadditive heritability that is not captured by LDSC when genetic interactions are indeed present in the generative model for a given complex trait. More importantly, iLDSC has a calibrated type I error rate and does not overestimate contributions of genetic interactions to trait variation in simulated data when only additive effects are present. While analyzing 25 complex traits in the UK Biobank and BioBank Japan, we illustrate that pairwise interactions are a source of ‘missing’ heritability captured by additive GWAS summary statistics—suggesting that phenotypic variation due to nonadditive genetic effects is more pervasive in human phenotypes than previously reported. Specifically, we find evidence of tagged genetic interaction effects contributing to heritability estimates in all of the 25 traits in the UK Biobank, and 23 of the 25 traits we analyzed in the BioBank Japan. We believe that iLDSC, with our development of a new cisinteraction score, represents a significant step towards resolving the true contribution of genetic interactions.
Results
Overview of the interactionLD score regression model
InteractionLD score regression (iLDSC) is a statistical framework for estimating heritability (i.e. the proportion of trait variance attributable to genetic variance). Here, we will give an overview of the iLDSC method and its corresponding software, as well as detail how its underlying model differs from that of LDSC (BulikSullivan et al., 2015b). We will assume that we are analyzing a GWAS dats set $\mathcal{D}=\{\mathbf{\mathbf{X}},\mathbf{\mathbf{y}}\}$ where $\mathbf{\mathbf{X}}$ is an $N\times J$ matrix of genotypes with $J$ denoting the number of SNPs (each of which is encoded as {0, 1, 2} copies of a reference allele at each locus $j$) and $\mathbf{\mathbf{y}}$ is an $N$dimensional vector of measurements of a quantitative trait. The iLDSC framework only requires summary statistics of individuallevel data: namely, marginal effect size estimates for each SNP $\widehat{\mathit{\bm{\beta}}}$ and a sample LD matrix $\mathbf{\mathbf{R}}$ (which can be provided via reference panel data).
We begin by considering the following generative linear model for complex traits
where ${b}_{0}$ is an intercept term; $\mathit{\bm{\beta}}=({\beta}_{1},\mathrm{\dots},{\beta}_{J})$ is a $J$dimensional vector containing the true additive effect sizes for an additional copy of the reference allele at each locus on $\mathbf{y}$; $\mathbf{W}$ is an $N\times M$ matrix of (pairwise) cisacting SNPbySNP statistical interactions between some subset of causal SNPs, where columns of this matrix are assumed to be the Hadamard (elementwise) product between genotypic vectors of the form ${\mathbf{\mathbf{x}}}_{j}\circ {\mathbf{\mathbf{x}}}_{k}$ for the $j$th and $k$th variants; $\mathit{\bm{\theta}}=({\theta}_{1},\mathrm{\dots},{\theta}_{M})$ is an $M$dimensional vector containing the interaction effect sizes; $\mathit{\bm{\epsilon}}$ is a normally distributed error term with mean zero and variance scaled according to the proportion of phenotypic variation not explained by genetic effects (BulikSullivan et al., 2015b), which we will refer to as the broadsense heritability of the trait denoted by ${H}^{2}$; and $\mathbf{\mathbf{I}}$ denotes an $N\times N$ identity matrix. For convenience, we will assume that the genotype matrix (columnwise) and the trait of interest have been meancentered and standardized (Strandén and Christensen, 2011; de Los Campos et al., 2013; Zhou et al., 2013). Lastly, we will let the intercept term b_{0} be a fixed parameter and we will assume that the effect sizes are each normally distributed with variances proportional to their individual contributions to trait heritability (Yang et al., 2010; Wu et al., 2011; Zhou et al., 2013; Crawford et al., 2017)
Effectively, we say that $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]={\phi}_{\beta}^{2}$ is the proportion of phenotypic variation contributed by additive SNP effects under the generative model, while $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]={\phi}_{\theta}^{2}$ makes up the proportion of phenotypic variation contributed by genetic interactions. While the appropriateness of treating genetic effects as random variables in analytical derivations has been questioned (de Los Campos et al., 2015), later, we will justify the theory presented here with simulation results showing that iLDSC accurately recovers nonadditive genetic variance in Equation 1 under a broad range of conditions.
There are two key takeaways from the generative model specified above. First, Equation 2 implies that the additive and nonadditive components in Equation 1 are orthogonal to each other. In other words, $\mathbb{E}[{\mathit{\bm{\beta}}}^{\u22ba}{\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=\mathbb{E}[{\mathit{\bm{\beta}}}^{\u22ba}]{\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{W}}\mathbb{E}[\mathit{\bm{\theta}}]=\mathrm{\U0001d7ce}$. This is important because it means that there is a unique partitioning of genetic variance when studying a trait of interest. The second key takeaway is that the genotype matrix $\mathbf{\mathbf{X}}$ and the matrix of genetic interactions $\mathbf{\mathbf{W}}$ themselves are correlated despite being linearly independent (see Materials and methods). This property stems from the fact that the pairwise interaction between two SNPs is encoded as the Hadamard product of two genotypic vectors in the form ${\mathbf{\mathbf{w}}}_{m}={\mathbf{\mathbf{x}}}_{j}\circ {\mathbf{\mathbf{x}}}_{k}$ (which is a nonlinear function of the genotypes).
A central objective in GWAS studies is to infer how much phenotypic variation can be explained by genetic effects. To achieve that objective, a key consideration involves incorporating the possibility of nonadditive sources of genetic variation to be explained by additive effect size estimates obtained from GWAS analyses (Hill et al., 2008). If we assume that the genotype and interaction matrices are correlated, then $\mathbf{X}$ and $\mathbf{\mathbf{W}}$ are not completely orthogonal (i.e. such that ${\mathbf{X}}^{\u22ba}\mathbf{W}\ne 0$) and the following relationship between the moment matrix $\mathbf{X}}^{\u22ba}\mathbf{y$, the observed marginal GWAS summary statistics $\hat{\beta}$, and the true coefficient values $\beta$ from the generative model in Equation 1 holds in expectation (see Materials and methods)
where $\mathbf{\mathbf{R}}$ is a sample estimate of the LD matrix, and $\mathbf{\mathbf{V}}$ represents a sample estimate of the correlation between the individuallevel genotypes $\mathbf{\mathbf{X}}$ and the span of genetic interactions between causal SNPs in $\mathbf{\mathbf{W}}$. Intuitively, the term $\mathbf{V}\theta$ can be interpreted as the subset of pairwise interaction effects that are tagged by the additive effect estimates from the GWAS study. Note that, when (i) nonadditive genetic effects do not contribute to the overall architecture of a trait (i.e. such that $\mathit{\theta}=\mathbf{0}$) or (ii) the genotype and interaction matrices $\mathbf{\mathbf{X}}$ and $\mathbf{\mathbf{W}}$ are uncorrelated, the equation above simplifies to a relationship between LD and summary statistics that is assumed in many GWAS studies and methods (Hormozdiari et al., 2014; Nakka et al., 2016; Zhu and Stephens, 2017; Zhang et al., 2018; Zhu and Stephens, 2018; Cheng et al., 2020; Demetci et al., 2021).
The goal of iLDSC is to increase estimates of genetic variance by accounting for sources of nonadditive genetic effects that can be explained by additive GWAS summary statistics. To do this, we extend the LD score regression framework and the corresponding LDSC software (BulikSullivan et al., 2015b). Here, according to Equation 3, we note that $\hat{\mathit{\beta}}\sim \mathcal{N}(\mathbf{R}\mathit{\beta}+\mathbf{V}\mathit{\theta},\lambda \mathbf{R})$ where $\lambda $ is a scale variance term due to uncontrolled confounding effects (Guan and Stephens, 2011; Song et al., 2022). Next, we condition on $\mathbf{\Theta}=(\mathit{\beta},\mathit{\theta})$ and take the expectation of chisquare statistics $\mathit{\chi}}^{2}=N\hat{\mathit{\beta}}{\hat{\mathit{\beta}}}^{\u22ba$ to yield
We define $\ell}_{j}=\sum _{k}{r}_{jk}^{2$ as the LD score for the additive effect of the $j$th variant (BulikSullivan et al., 2015b), and $f}_{j}=\sum _{m}{v}_{jm}^{2$ represents the ‘cisinteraction’ LD score which encodes the pairwise interaction between the $j$th variant and all other variants within a genomic window that is a prespecified number of SNPs wide (Crawford et al., 2017), respectively. By considering only the diagonal elements of LD matrix in the first term, similar to the original LDSC approach (BulikSullivan et al., 2015b; Song et al., 2022), we get the following simplified regression model
where ${\mathit{\chi}}^{2}=({\chi}_{1}^{2},\dots ,{\chi}_{J}^{2})$ is a $J$dimensional vector of chisquare summary statistics, and $\mathit{\ell}=({\ell}_{1},\dots ,{\ell}_{J})$ and $\mathit{f}=({f}_{1},\dots ,{f}_{J})$ are $J$dimensional vectors of additive and cisinteraction LD scores, respectively. Furthermore, we define the variance components $\tau =N{\phi}_{\beta}^{2}/J$ and $\vartheta =N{\phi}_{\theta}^{2}/M$ as the additive and nonadditive regression coefficients of the model, and 1 is the intercept meant to model the bias factor due to uncontrolled confounding effects (e.g. cryptic relatedness structure). In practice, we efficiently compute the cisinteraction LD scores by considering only a subset of interactions between each $j$th focal SNP and SNPs within a cisproximal window around the $j$th SNP. In our validation studies and applications, we base the width of this window on the observation that LD decays outside of a window of 1 centimorgan (cM); therefore, SNPs outside the 1 cM window centered on the $j$th SNP will not significantly contribute to its LD scores. Note that the width of this window can be relaxed in the iLDSC software when appropriate. We fit the iLDSC model using weighted least squares to estimate regression parameters and derive pvalues for identifying traits that have significant statistical evidence of tagged cisinteraction effects by testing the null hypothesis ${H}_{0}:\vartheta =0$. Importantly, under the null model of a trait being generated by only additive effects, the iLDSC model in Equation 5 reduces to an infinitesimal model (Fisher, 1999) or, in the case some variants have no effect on the trait, a polygenic model.
Lastly, we want to note the empirical observation that the additive ($\mathbf{\ell}$) and interaction $(\mathit{\bm{f}}$) LD scores are lowly correlated. This is important because it indicates that the presence of cisinteraction LD scores in the model specified in Equation 5 has littletono influence over the estimate for the additive coefficient $\tau $. Instead, the inclusion of $\mathit{\bm{f}}$ creates a multivariate model that can identify the proportion of variance explained by both additive and nonadditive effects in summary statistics. In other words, we can interpret $\widehat{\vartheta}$ as an estimate of the phenotypic variation explained by tagged cisacting interaction effects. The concept of additive genetic effects partially explaining nonadditive variation has also described in various studies from quantitative genetics (Hill et al., 2008; Hivert et al., 2021; MäkiTanila and Hill, 2014). Under HardyWeinberg equilibrium, it can be shown that the additive variance explained by $J$ SNPs takes on the following form (Materials and methods) (Falconer and Mackay, 1983)
The expression for the additive variance ${\sigma}_{A}^{2}$ in Equation 6 is important because it represents the theoretical upper bound on the proportion of total phenotypic variance that can be recovered from GWAS summary statistics using the iLDSC framework. As a result, we use the sum of coefficient estimates $\widehat{\tau}+\widehat{\vartheta}\le {\sigma}_{A}^{2}$ to construct iLDSC heritability estimates. A full derivation of the cisinteraction regression framework and details about its corresponding implementation in our software iLDSC can be found in Materials and Methods.
Detection of tagged pairwise interaction effects using iLDSC in simulations
We illustrate the power of iLDSC across different genetic trait architectures via extensive simulation studies (Materials and methods). We generate synthetic phenotypes using real genomewide genotype data from individuals of selfidentified European ancestry in the UK Biobank. To do so, we first assume that traits have a polygenic architecture where all SNPs have a nonzero additive effect. Next, we randomly select a set of causal cisinteraction variants and divide them into two interacting groups (Materials and methods). One may interpret the SNPs in group #1 as being the ‘hubs’ in an interaction map (Crawford et al., 2017), whereas SNPs in group #2 are selected to be variants within some kilobase (kb) window around each SNP in group #1. We assume a wide range of simulation scenarios by varying the following parameters:
heritability: H^{2} = 0.3 and 0.6;
proportion of phenotypic variation that is generated by additive effects: $\rho =$ 0.5, 0.8, and 1;
percentage of SNPs selected to be in group #1: 1%, 5%, and 10%;
genomic window used to assign SNPs to group #2: ± 10 and ± 100 kb.
We also varied the correlation between SNP effect size and minor allele frequency (MAF; as discussed in Schoech et al., 2019). All results presented in this section are based on 100 different simulated phenotypes for each parameter combination.
Figure 1 demonstrates that iLDSC robustly detects significant tagged nonadditive genetic variance, regardless of the total number of causal interactions genomewide. Instead, the power of iLDSC depends on the proportion of phenotypic variation that is generated by additive versus interaction effects (ρ), and its power tends to scale with the window size used to compute the cisinteraction LD scores (see Materials and methods). iLDSC shows a similar performance for detecting tagged cisinteraction effects when the effect sizes of causal SNPs depend on their minor allele frequency and when we varied the number of SNPs assigned to be in group #2 within 10 kb and 100 kb windows, respectively (Figure 1—figure supplements 1–5).
Importantly, iLDSC does not falsely identify putative nonadditive genetic effects in GWAS summary statistics when the synthetic phenotype was generated by only additive effects ($\rho =1$). Figure 2 illustrates the performance of iLDSC under the null hypothesis ${H}_{0}:\vartheta =0$, with the type I error rates for different estimation window sizes of the cisinteraction LD scores highlighted in panel A. Here, we also show that, when no genetic interaction effects are present, iLDSC unbiasedly estimates the cisinteraction coefficient in the regression model to be $\widehat{\vartheta}=0$ (Figure 2B), robustly estimates the heritability (Figure 2C), and provides wellcalibrated pvalues when assessed over many traits (Figure 2D). This behavior is consistent across different MAFdependent effect size distributions, and pvalue calibration is not sensitive to misspecification of the estimation windows used to generate the cisinteraction LD scores (Figure 2—figure supplements 1–2).
One of the innovations that iLDSC offers over the traditional LDSC framework is increased heritability estimates after the identification of nonadditive genetic effects that are tagged by GWAS summary statistics. Here, we applied both methods to the same set of simulations in order to understand how LDSC behaves for traits generated with cisinteraction effects. Figure 3 depicts boxplots of the heritability estimates for each approach and shows that, across an array of different synthetic phenotype architectures, LDSC captures less of phenotypic variance explained by all genetic effects. It is important to note that iLDSC can yield upwardly biased heritability estimates when the cisinteraction scores are computed over genomic window sizes that are too small; however, these estimates become more accurate for larger window size choices (Figure 3—figure supplement 1). In contrast to LDSC, which aims to capture phenotypic variance attributable to the additive effects of genotyped SNPs, iLDSC accurately partitions genetic effects into additive versus cisinteracting components, which in turn generally leads the ability of iLDSC to capture more genetic variance. The mean absolute error between the true generative heritability and heritability estimates produced by iLDSC and LDSC are shown in Supplementary files 1 and 2, respectively. Generally, the error in heritability estimates is higher for LDSC than it is for iLDSC across each of the scenarios that we consider.
Next, we perform an additional set of simulations where we explore other common generative models for complex trait architecture that involve nonadditive genetic effects. Specifically, we compare heritability estimates from LDSC and iLDSC in the presence of additive effects, cisacting interactions, and a third source of genetic variance stemming from either genebyenvironment (G×E) or or genebyancestry (G×Ancestry) effects. Details on how these components were generated can be found in Materials and Methods. In general, iLDSC underestimates overall heritability when additive effects and cisacting interactions are present alongside G×E (Figure 3—figure supplement 2) and/or G×Ancestry effects when PCs are included as covariates (Figure 3—figure supplement 3). Notably, when PCs are not included to correct for residual stratification, both LDSC and iLDSC can yield unbounded heritability estimates greater than 1 (Figure 3—figure supplement 4). Also interestingly, when we omit cisinteractions from the generative model (i.e. the genetic architecture of simulated traits is only made up of additive and G×E or G×Ancestry effects), iLDSC will still estimate a nonzero genetic variance component with the cisinteraction LD scores (Figure 3—figure supplements 5–7). Collectively, these results empirically show the important point that cisinteraction scores are not enough to recover missing genetic variation for all types of trait architectures; however, they are helpful in recovering phenotypic variation explained by statistical interaction effects. Recall that the linear relationship between (expected) ${\chi}^{2}$ test statistics and LD scores proposed by the LDSC framework holds when complex traits are generated under the polygenic model where all causal variants have the same expected contribution to phenotypic variation. When cisinteractions affect genetic architecture (e.g. in our earlier simulations in Figure 3), these assumptions are violated in LDSC, but the inclusion of the additional nonlinear scores in iLDSC help recover the relationship between the expectation of ${\chi}^{2}$ test statistics and LD.
As a further demonstration of how iLDSC performs when assumptions of the original LD score model are violated, we also generated synthetic phenotypes with sparse architectures using the spikeandslab model (Zhou et al., 2013). Here, traits were simulated with solely additive effects, but this time only variants with the top or bottom $\{1,5,10,25,50,100\}$ percentile of LD scores were given nonzero effects (see Materials and methods). Breaking the relationship assumed under the LDSC framework between LD scores and chisquared statistics (i.e. that they are generally positively correlated) led to unbounded estimates of heritability in all but the (polygenic) scenario when 100% of SNPs contributed to the phenotypic variation (Figure 3—figure supplement 8).
Finally, we performed a set of polygenic simulations to assess if iLDSC estimates of nonadditive genetic variance could be spuriously inflated due to either (i) unobserved additive effects (see, for example, Hemani et al., 2014), (ii) unobserved SNPs that are involved in genetic interactions, or by (iii) nonzero correlation between the additive and interaction effect sizes in the generative model (i.e. breaking the independence assumption in Equation 2). In the first setting, we observed that, across a range of both minor allele frequencies and effect sizes, the omission of causal haplotypes had a negligible effect on the estimated value of the coefficients in iLDSC (Figure 3—figure supplement 9). We hypothesize this is due to the fact that the simulations were done for polygenic architectures where all SNPs have at least an additive effect. As a result, not observing a small subset of SNPs does not hinder the ability of iLDSC to estimate genetic variance because the effect size of each SNP is small. If these simulations were conducted for sparse architectures, we would have likely seen a greater impact on iLDSC; although, we have already shown the LD score regression framework to be uncalibrated for traits with sparse genetic architectures (again see Figure 3—figure supplement 8). In the second setting, we observed that the iLDSC framework protects against the false discovery of nonadditive genetic effects and underestimates the variance component $\vartheta $ when causal variants involved in pairwise interactions were unobserved (Figure 3—figure supplements 10 and 11). As a direct comparison, estimates of the additive variance component $\tau $ in iLDSC were not affected by the unobserved interacting variants. Lastly, in the third setting, we observed that the mean estimate of the genetic variance in both LDSC and iLDSC had a slight upward bias as the correlation between additive and interaction effect sizes in the generative model increased; however, the median of these bias estimates was still near zero across all simulated scenarios and their corresponding replicates (Figure 3—figure supplements 12 and 13).
Application of iLDSC to the UK Biobank and BioBank Japan
To assess whether pairwise interaction genetic effects are significantly affecting estimates of heritability in empirical biobank data, we applied iLDSC to 25 continuous quantitative traits from the UK Biobank and BioBank Japan (Supplementary file 3). Protocols for computing GWAS summary statistics for the UK Biobank are described in the Materials and methods; while precomputed summary statistics for BioBank Japan were downloaded directly from the consortium website (https://pheweb.jp/downloads). We release the cisacting SNPbySNP interaction LD scores used in our analyses on the iLDSC GitHub repository from two reference groups in the 1000 Genomes: 489 individuals from the European superpopulation (EUR) and 504 individuals from the East Asian (EAS) superpopulation (see also Supplementary files 4 and 5).
In each of the 25 traits, we analyzed in the UK Biobank, we detected significant proportions of estimated genetic variation stemming from tagged pairwise cisinteractions (Table 1). This includes many canonical traits of interest in heritability analyses: height, cholesterol levels, urate levels, and both systolic and diastolic blood pressure. Our findings in Table 1 are supported by multiple published studies identifying evidence of nonadditive effects playing a role in the architectures of different traits of interest. For example, Li et al., 2020 found evidence for genetic interactions that contributed to the pathogenesis of coronary artery disease. It was also recently shown that nonadditive genetic effects plays a significant role in body mass index (Song et al., 2022). Generally, we find that the traditional LDSC produces lower estimates of trait heritability because it does not consider the additional sources of genetic signal that iLDSC does (Table 1). In BioBank Japan, 23 of the 25 traits analyzed had a significant nonlinear component detected by iLDSC — with HDL and triglyceride levels being the only exceptions.
For each of the 25 traits that we analyzed, we found that the iLDSC heritability estimates are significantly correlated with corresponding estimates from LDSC in both the UK Biobank (${r}^{2}=0.988$, $P=5.936\times {10}^{24}$) and BioBank Japan (${r}^{2}=0.849$, $P=6.061\times {10}^{11}$) as shown in Figure 4A. Additionally, we found that the heritability estimates for the same traits between the two biobanks are highly correlated according to both LDSC (${r}^{2}=0.848$, $P=7.166\times {10}^{11}$) and iLDSC (${r}^{2}=0.666$, $P=6.551\times {10}^{7}$) analyses as shown in Figure 4B. After comparing the iLDSC heritability estimates to LDSC, we then assessed whether there was significant difference in the amount of phenotypic variation explained by the nonadditive genetic effect component in the GWAS summary statistics derived from the the UK Biobank and BioBank Japan (i.e. comparing the estimates of $\vartheta $; see Figure 4—figure supplement 1A). We show that, while heterogeneous between traits, the phenotypic variation explained by genetic interactions is relatively of the same magnitude for both biobanks (${r}^{2}=0.372$, $P=0.0119$). Notably, the trait with the most significant evidence of tagged cisinteraction effects in GWAS summary statistics is height which is known to have a highly polygenic architecture.
The intercepts estimated by LDSC and iLDSC are also highly correlated in both the UK Biobank and the BioBank Japan (Figure 4—figure supplement 1B). Recall that these intercept estimates represent the confounding factor due to uncontrolled effects. For LDSC, this does include phenotypic variation that is due to unaccounted for pairwise statistical genetic interactions. The iLDSC intercept estimates tend to be correlated with, but are generally different than, those computed with LDSC — empirically indicating that nonadditive genetic variation is partitioned away and is missed when using the standard LD score alone. This result shows similar patterns in both the UK Biobank (${r}^{2}=0.888$, $P=1.962\times {10}^{12}$) and BioBank Japan (${r}^{2}=0.813$, $P=7.814\times {10}^{10}$).
Lastly, we performed an additional analysis in the UK Biobank where the cisinteraction scores are included as an annotation alongside 97 other functional categories in the stratifiedLD score regression framework and its software sLDSC (Gazal et al., 2017; Materials and methods). Here, sLDSC heritability estimates still showed an increase with the interaction scores versus when the publicly available functional categories were analyzed alone, but albeit at a much smaller magnitude (Table 2). The contributions from the pairwise interaction component to the overall estimate of genetic variance ranged from 0.005 for MCHC ($P=0.373$) to 0.055 for HDL ($P=0.575$; Figure 4C and D). Furthermore, in this analysis, the estimates of the nonadditive components were no longer statistically significant for any of the traits in the UK Biobank (Table 2). Despite this, these results highlight the ability of the iLDSC framework to identify sources of ‘missing’ phenotypic variance explained in heritability estimation. Importantly, moving forward, we suggest using the cisinteraction scores with additional annotations whenever they are available as it provides more conservative estimates of the role of nonadditive effects on trait architecture.
Discussion
In this paper, we present iLDSC, an extension of the LD score regression framework which aims to recover missing heritability from GWAS summary statistics by incorporating an additional score that measures the nonadditive genetic variation that is tagged by genotyped SNPs. Here, we demonstrate how iLDSC builds upon the original LDSC model through the development of new ‘cisinteraction’ LD scores which help to investigate signals of cisacting SNPbySNP interactions (Figure 1 and Figure 1—figure supplements 1–5). Through extensive simulations, we show that iLDSC is wellcalibrated under the null model when polygenic traits are generated only by additive effects (Figure 2 and Figure 2—figure supplements 1–2), we highlight that iLDSC provides greater heritability estimates over LDSC when traits are indeed generated with cisacting SNPbySNP interaction effects (Figure 3 and Figure 3—figure supplement 1, and Supplementary files 1 and 2), and we tested the robustness of iLDSC on phenotypes where assumptions of the original LD score model are violated (Figure 3—figure supplements 2–13). Finally, in real data, we show examples of many traits with estimated GWAS summary statistics that tag cisinteraction effects in the UK Biobank and BioBank Japan (Figure 4 and Figure 4—figure supplement 1, Tables 1 and 2, and Supplementary files 35). We have made iLDSC a publicly available command line tool that requires minimal updates to the computing environment used to run the original implementation of LD score regression. In addition, we provide precomputed cisinteraction LD scores calculated from the European (EUR) and East Asian (EAS) reference populations in the 1000 Genomes phase 3 data (see Data and Software Availability under Materials and Methods).
The current implementation of the iLDSC framework offers many directions for future development and applications. First, an area of future work would be to explore how the relationship between cisinteraction LD scores and interaction effect sizes from the generative model of complex traits might bias heritability estimates provided by iLDSC (e.g., similar to the relationship we explored between the standard LD scores and linear effect sizes in Figure 3—figure supplement 8). Second, as we showed with our simulation studies (Figure 3—figure supplements 2–8), the cisinteraction LD scores that we propose are not always enough to recover explainable nonadditive genetic effects for all types of trait architectures. While we focus on pairwise cisacting SNPbySNP statistical interactions in this work, the theoretical concepts underlying iLDSC can easily be adapted to other types of interactions as well. Third, in our analysis of the UK Biobank and BioBank Japan, we showed that the inclusion of additional categories via frameworks such as stratified LD score regression (Finucane et al., 2015) can be used to provide more refined heritability estimates from GWAS summary statistics while accounting for linkage (see results in Table 1 versus Table 2). A key part of our future work is to continue to explore whether considering functional annotation groups would also improve our ability to identify tagged nonadditive genetic effects. Lastly, we have only focused on analyzing one phenotype at a time in this study. However, many previous studies have extensively shown that modeling multiple phenotypes can often dramatically increase power (Runcie et al., 2020; Stamp et al., 2022). Therefore, it would be interesting to extend the iLDSC framework to multiple traits to study nonlinear genetic correlations in the same way that LDSC was recently extended to uncover additive genetic correlation maps across traits (Naqvi et al., 2021).
Materials and methods
Generative statistical model for complex traits
Request a detailed protocolOur goal in this study is to reanalyze summary statistics from genomewide association studies (GWAS) and estimate heritability while accounting for both additive genetic associations and tagged interaction effects. We begin by assuming the following generative linear model for complex traits which can be seen as an extended view of Equation 1 in the main text
where $\mathbf{\mathbf{y}}$ denotes an $N$dimensional vector of phenotypic states for a quantitative trait of interest measured in $N$ individuals; ${b}_{0}$ is an intercept term; $\mathbf{\mathbf{X}}$ is an $N\times J$ matrix of genotypes, with $J$ denoting the number of single nucleotide polymorphism (SNPs) encoded as $\{0,1,2\}$ copies of a reference allele at each locus; $\mathit{\bm{\beta}}=({\beta}_{1},\mathrm{\dots},{\beta}_{J})$ is a $J$dimensional vector containing the true additive effect sizes for an additional copy of the reference allele at each locus on $\mathbf{y}$; $X}_{D$ is an $N\times J$ matrix that represents the dominance for each genotype encoded as $\{0,1,1\}$ with corresponding effect sizes $\omega$; $\mathbf{W}$ is an $N\times M$ matrix of genetic interactions; $\mathit{\bm{\theta}}=({\theta}_{1},\mathrm{\dots},{\theta}_{M})$ is an $M$dimensional vector containing the interaction effect sizes; $\mathit{\bm{\epsilon}}$ is a normally distributed error term with mean zero and variance scaled according to the proportion of phenotypic variation not explained by the broadsense heritability of the trait, denoted by ${H}^{2}$; and $\mathbf{\mathbf{I}}$ denotes an $N\times N$ identity matrix. Note that the encoding for dominance in ${\mathbf{\mathbf{X}}}_{D}$ was chosen because it imposes orthogonality with the genotype encoding in $\mathbf{\mathbf{X}}$ (Purcell et al., 2007; Vitezica et al., 2017; Palmer et al., 2023).
For convenience, we will assume that the genotype matrix (columnwise), the dominance matrix (also columnwise), and trait of interest have all been standardized (Strandén and Christensen, 2011; de Los Campos et al., 2013; Zhou et al., 2013). Furthermore, while the matrix $\mathbf{\mathbf{W}}$ could encode any source of nonadditive genetic interactions (e.g. genebyenvironmental effects) in theory, we limit our focus in this study to trait architectures that have been generated with contributions stemming from cisacting statistical SNPbySNP (or pairwise) interactions. To that end, we assume that the columns of $\mathbf{\mathbf{W}}$ are the Hadamard (elementwise) product between genotypic vectors of the form ${\mathbf{\mathbf{x}}}_{j}\circ {\mathbf{\mathbf{x}}}_{k}$ for the $j$th and $k$th variants. We also want to point out that the generative formulation of Equation 7 can also be easily extended to accommodate other fixed effects (e.g. age, sex, or genotype principal components), as well as other random effects terms that can be used to account for sample nonindependence due to other environmental factors.
As a final set of assumptions, we will let the intercept term ${b}_{0}$ be a fixed parameter while allowing the other coefficients to follow independent Gaussian distributions with variances proportional to their individual contributions to the trait heritability (Yang et al., 2010; Wu et al., 2011; Zhou et al., 2013; Jiang and Reif, 2015; Crawford et al., 2017)
for $j=1,\mathrm{\dots},J$ and $m=1,\mathrm{\dots},M$. The broadsense heritability of the trait is defined as ${H}^{2}={\phi}_{\beta}^{2}+{\phi}_{\omega}^{2}+{\phi}_{\theta}^{2}$. Under the generative model in Equation 7, we then say that $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]={\phi}_{\beta}^{2}$ is the proportion of phenotypic variation contributed by additive SNP effects, $\mathbb{V}[{\mathbf{\mathbf{X}}}_{D}\mathit{\bm{\omega}}]={\phi}_{\omega}^{2}$ is the proportion of phenotypic variation contributed by dominance effects, and the set of interactions involving some subset of causal SNPs contribute the remaining proportion to the heritability $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]={\phi}_{\theta}^{2}$. As we mentioned in the main text, we recognize that the appropriateness of treating genetic effects as random variables in analytical derivations has been questioned (de Los Campos et al., 2015), but our simulation studies show that iLDSC accurately recovers nonadditive genetic variance in Equation 7 under a broad range of conditions.
Orthogonality between additive and nonadditive genetic effects
Request a detailed protocolAssuming that the effect sizes $\{\mathit{\bm{\beta}},\mathit{\bm{\omega}},\mathit{\bm{\theta}}\}$ in Equation 8 follow independent and zero mean Gaussian distributions leads to orthogonality between the additive and nonadditive components in Equation 7. Since the genotypes $\mathbf{\mathbf{X}}$ and the dominance values ${\mathbf{\mathbf{X}}}_{D}$ are fixed orthogonal matrices, it is straightforward to show that $\text{Cov}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}},{\mathbf{\mathbf{X}}}_{D}\mathit{\bm{\omega}}]=0$ (Vitezica et al., 2017; Palmer et al., 2023). The same relationship can be shown for the additive and the pairwise interaction genetic effects where
with ${\mathbf{\mathbf{x}}}_{j}$ and ${\mathbf{\mathbf{w}}}_{m}$ denoting the $j$th and $m$th column of the individuallevel genotype matrix $\mathbf{\mathbf{X}}$ and the interaction matrix $\mathbf{\mathbf{W}}$, respectively. Note that a similar derivation to Equation 9 can also be done for the dominance and pairwise genetic interaction effects. This concept of orthogonality is important because we want to preserve a unique partitioning of genetic variance when modeling a trait of interest.
Genotypes and their interactions are correlated despite being linearly independent
The design matrices $\mathbf{\mathbf{X}}$ and $\mathbf{\mathbf{W}}$ in Equation 7 are not linearly dependent because the pairwise interactions between two SNPs are encoded as the Hadamard product of two genotypic vectors in the form ${\mathbf{\mathbf{x}}}_{j}\circ {\mathbf{\mathbf{x}}}_{k}$ (which is a nonlinear function). Linear dependence would have implied that one could find a transformation between a SNP and an interaction term in the form ${\mathbf{\mathbf{w}}}_{m}=c\times {\mathbf{\mathbf{x}}}_{j}$ for some constant $c$. However, despite their linear independence, $\mathbf{\mathbf{X}}$ and $\mathbf{\mathbf{W}}$ are themselves not orthogonal and still have a nonzero correlation. This implies that the inner product between genotypes and their interactions is nonzero ${\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{W}}\ne \mathrm{\U0001d7ce}$. To see this, we focus on a focal SNP ${\mathbf{\mathbf{x}}}_{j}$ and consider three different types of interactions:
Scenario I: Interaction between a focal SNP with itself (${\mathbf{\mathbf{x}}}_{j}\circ {\mathbf{\mathbf{x}}}_{j}$).
Scenario II: Interaction between a focal SNP with a different SNP (${\mathbf{\mathbf{x}}}_{j}\circ {\mathbf{\mathbf{x}}}_{k}$).
Scenario III: Interaction between a focal SNP with a pair of different SNPs (${\mathbf{\mathbf{x}}}_{k}\circ {\mathbf{\mathbf{x}}}_{l}$).
The following derivations rely on the fact that: (1) we assume that genotypes have been meancentered and scaled to have unit variance, and (2) under HardyWeinberg equilibrium, SNPs marginally follow a binomial distribution ${\mathbf{\mathbf{x}}}_{j}\sim \text{Bin}(2,p)$ where $p$ represents the minor allele frequency (MAF) (Wray et al., 2007; Lippert et al., 2013).
Scenario I
Request a detailed protocolThe covariance between a focal SNP and an interaction with itself is $\text{Cov}[{\mathbf{\mathbf{x}}}_{j},{\mathbf{\mathbf{x}}}_{j}{\mathbf{\mathbf{x}}}_{j}]=\mathbb{E}[{\mathbf{\mathbf{x}}}_{j}^{3}]\mathbb{E}[{\mathbf{\mathbf{x}}}_{j}]\mathbb{E}[{\mathbf{\mathbf{x}}}_{j}^{2}]$. With meancentered SNPs, this is proportional to $\mathbb{E}[{\mathbf{\mathbf{x}}}_{j}^{3}]=(qp)/\sqrt{2pq}$ which is the skewness of the binomial distribution where, again, $p=$ MAF and $q=$ 1MAF of the $j$th SNP.
Scenario II
Request a detailed protocolAssume that we have two SNPs, ${\mathbf{\mathbf{x}}}_{j}\sim \text{Bin}(2,{p}_{j})$ and ${\mathbf{\mathbf{x}}}_{k}\sim \text{Bin}(2,{p}_{k})$ where ${p}_{j}$ and ${p}_{k}$ represent their respective minor allele frequencies. We want to compute the correlation between ${\mathbf{\mathbf{x}}}_{j}$ and the interaction ${\mathbf{\mathbf{x}}}_{j}{\mathbf{\mathbf{x}}}_{k}$ where $\text{Cov}[{\mathbf{\mathbf{x}}}_{j},{\mathbf{\mathbf{x}}}_{j}{\mathbf{\mathbf{x}}}_{k}]=\mathbb{E}[{\mathbf{\mathbf{x}}}_{j}^{2}{\mathbf{\mathbf{x}}}_{k}]\mathbb{E}[{\mathbf{\mathbf{x}}}_{j}]\mathbb{E}[{\mathbf{\mathbf{x}}}_{j}{\mathbf{\mathbf{x}}}_{k}]$. Again, with the meancentered assumption, the covariance is proportional to the expectation $\mathbb{E}[{\mathbf{\mathbf{x}}}_{j}^{2}{\mathbf{\mathbf{x}}}_{k}]$. Here, with SNPs taking on values $\{0,1,2\}$, the joint distribution between ${\mathbf{\mathbf{x}}}_{j}^{2}$ and ${\mathbf{\mathbf{x}}}_{k}$ can be written out as the following Kang and Jung, 2001:
${\mathbf{\mathbf{x}}}_{j}^{2}=0$  ${\mathbf{\mathbf{x}}}_{j}^{2}=1$  ${\mathbf{\mathbf{x}}}_{j}^{2}=4$  

${\mathbf{\mathbf{x}}}_{k}=0$  ${u}_{jk}^{2}$  $2{u}_{jk}(1{p}_{k}{u}_{jk})$  ${(1{p}_{k}{u}_{jk})}^{2}$ 
${\mathbf{\mathbf{x}}}_{k}=1$  $2{u}_{jk}(1{p}_{j}{u}_{jk})$  $\begin{array}{l}2{u}_{jk}({u}_{jk}+{p}_{j}+{p}_{k}1)+\\ 2(1{p}_{j}{u}_{jk})(1{p}_{k}{u}_{jk})\end{array}$  $2({u}_{jk}+{p}_{j}+{p}_{k}1)(1{p}_{k}{u}_{jk})$ 
${\mathbf{\mathbf{x}}}_{k}=2$  ${(1{p}_{j}{u}_{jk})}^{2}$  $2({u}_{jk}+{p}_{j}+{p}_{k}1)(1{p}_{j}{u}_{jk})$  ${({u}_{jk}+{p}_{j}+{p}_{k}1)}^{2}$ 
where ${u}_{jk}=(1{p}_{j})(1{p}_{k})+{r}_{jk}\sqrt{{p}_{j}{p}_{k}(1{p}_{j})(1{p}_{k})}$ and ${r}_{jk}$ is the Pearson correlation or linkage disequilibrium (LD) between the $j$th and $k$th SNPs.
Scenario III
Request a detailed protocolThe covariance between a focal SNP and an interaction with a pair of different SNPs $\text{Cov}[{\mathbf{\mathbf{x}}}_{j},{\mathbf{\mathbf{x}}}_{k}{\mathbf{\mathbf{x}}}_{l}]$ will be nonzero if the $j$th SNP is correlated with either variant (i.e., ${r}_{jk}\ne 0$ or ${r}_{jl}\ne 0$).
Traditional estimation of additive GWAS summary statistics
Request a detailed protocolAs previously mentioned, the key to this work is that SNPlevel GWAS summary statistics can also tag nonadditive genetic effects when there is a nonzero correlation between individuallevel genotypes and their interactions (as defined in Equation 7). Throughout the rest of this section, we will use ${\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{X}}/N$ to denote the LD or pairwise correlation matrix between SNPs. We will then let $\mathbf{\mathbf{R}}$ represent an LD matrix empirically estimated from external data (e.g. directly from GWAS study data, or using a pairwise LD map from a population that is representative of the samples analyzed in the GWAS study). The important property here is the following
where the term ${r}_{jk}$ is again defined as the Pearson correlation coefficient between the $j$th and $k$th SNPs, respectively.
In traditional GWAS studies, summary statistics of the true additive effects $\mathit{\bm{\beta}}={({\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{X}})}^{1}{\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{y}}$ in Equation 7 are typically derived by computing a marginal least squares estimate with the observed data
There are two key identities that may be taken from Equation 11. The first uses Equation 10 and is the approximate relationship (in expectation) between the moment matrix ${\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{y}}$ and the linear effect size estimates $\widehat{\mathit{\bm{\beta}}}$:
The second key point combines Equations 10 and 12 to describe the asymptotic relationship between the observed marginal GWAS summary statistics $\widehat{\mathit{\bm{\beta}}}$ and the joint coefficient values $\mathit{\bm{\beta}}$ where (in expectation)
After some algebra, the above mirrors a highdimensional regression model (in expectation) where $\widehat{\mathit{\bm{\beta}}}=\mathbf{\mathbf{R}}\mathit{\bm{\beta}}$ with the estimated summary statistics as the response variables and the empirically estimated LD matrix acting as the design matrix (Hormozdiari et al., 2014; Hormozdiari et al., 2016; Zhang et al., 2018; Cheng et al., 2020; Demetci et al., 2021). Theoretically, the resulting coefficients output from this highdimensional model are the desired true effect size estimates used to generate the phenotype of interest.
Additive GWAS summary statistics with tagged interaction effects
Request a detailed protocolWhen interactions contribute to the architecture of complex traits (i.e. $\mathit{\bm{\theta}}\ne \mathrm{\U0001d7ce}$), the marginal GWAS summary statistics derived using least squares in Equation 11 will also explain nonadditive variation when there is a nonzero correlation between genotypes and their interactions. To see this, we use the concept of ‘omitted variable bias’ (Barreto and Howland, 2005) where the fitted model aims to estimate the true additive coefficients $\mathit{\bm{\beta}}$ but does not account for contributions from the nonadditive components which also contribute to trait architecture. In this case, we get the following
Since we assume that the genotypes are orthogonal to both the dominance effects in Equation 7, we know that ${\mathbf{\mathbf{X}}}^{\u22ba}{\mathbf{\mathbf{X}}}_{D}=\mathrm{\U0001d7ce}$. This simplifies the above to be the following
where the matrix ${\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{W}}$(which we showed to be nonzero) can be interpreted as the sample correlation between individuallevel genotypes and the cisinteractions between causal SNPs. By taking the expectation using Equations 10 and 12, we get the following alternative (approximate) relationship between the observed marginal GWAS summary statistics $\widehat{\mathit{\bm{\beta}}}$ and the true coefficient values $\mathit{\bm{\beta}}$
which results from our initial assumption that the residuals are normally distributed with mean zero $\mathbb{E}[\mathit{\bm{\epsilon}}]=\mathrm{\U0001d7ce}$ in Equation 7. Here, we define $\mathbf{\mathbf{V}}$ to represent a sample estimate of the correlation between the individuallevel genotypes and the nonadditive genetic interaction matrix such that $\mathbb{E}[{\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{W}}]\approx N\mathbf{\mathbf{V}}$. Similar to the LD matrix $\mathbf{\mathbf{R}}$, the correlation matrix $\mathbf{\mathbf{V}}$ is also assumed to be computed from reference panel data. Intuitively, when $\mathit{\bm{\theta}}\ne \mathrm{\U0001d7ce}$ there is additional phenotypic variation contributed by pairwise interactions that can be explained by GWAS effect size estimates. Moreover, when $\mathbf{\mathbf{V}}\mathit{\bm{\theta}}=\mathrm{\U0001d7ce}$, then the relationship in Equation 16 converges onto the conventional asymptotic assumption (in expectation) between GWAS summary statistics and the true additive coefficients in Equation 13; Hormozdiari et al., 2014; Hormozdiari et al., 2016; Zhang et al., 2018; Cheng et al., 2020; Demetci et al., 2021.
Connection to quantitative genetics theory
Request a detailed protocolThe concept of additive genetic effects partially explaining nonadditive variation has also described in classical quantitative genetics (Hill et al., 2008; Hivert et al., 2021; MäkiTanila and Hill, 2014). Consider an individual genotyped at $J$ loci each with major and minor alleles A and B, respectively. Let ${p}_{j}$ be the allele frequency of A at the $j$th locus, ${a}_{j}$ denote the additive effect, and ${[aa]}_{jk}$ be the additivebyadditive (pairwise) interaction effect between loci $j$ and $k$, and ${[aaa]}_{jkl}$ represent a third order interaction between loci $j$, $k$, and $l$. For simplicity in presentation, assume that dominance only makes a small contribution to the genetic variance (Palmer et al., 2023; Pazokitoroudi et al., 2021; Zhu et al., 2015). The population mean is given as the following
We follow the assumption that the genetic variation in human complex traits can predominately be explained by additive effects, with the remainder variation being mostly explained by additivebyadditive effects (Weinreich et al., 2018; Jiang and Reif, 2015; Fisher, 1919; Lynch and Walsh, 1998). As a result, we will ignore the higher order interaction terms in Equation 17. Under HardyWeinberg equilibrium, we can find the average effect by taking the first derivative of the population mean with respect to the frequency of the increasing allele (MäkiTanila and Hill, 2014; Hivert et al., 2021). For the $j$th SNP, the average effect (including terms up to secondorder interaction) is given by the following
which notably contains both the additive effect and a summation of additivebyadditive interactions between pairs of loci. The additive genetic variance for the $j$th SNP takes on the following form
which is the product of the square of the average effect in Equation 18 and the heterozygosity at $j$th locus $\mathbb{V}[{\mathbf{\mathbf{x}}}_{j}]=2{p}_{j}(1{p}_{j})$ (again assuming that SNPs marginally follow a binomial distribution). The total additive variance is then obtained by summing over the $J$ loci such that ${\sigma}_{A}^{2}={\sum}_{j}{\sigma}_{A}^{2}(j)$ (Falconer and Mackay, 1983).
We can derive a parallel construction for additive genetic variance using the generative random effect model presented in Equation 7; Hivert et al., 2021. Here, we will leverage that with genotype data taken for $N$ individuals, ${\sum}_{i}{x}_{ij}/N=2{p}_{j}$. Ignoring the assumed small contributions from dominance effects, the population mean for a quantitative trait $\mathbf{\mathbf{y}}$ can be written as the following
To find the average effect for the $j$th locus, we this time take the first derivative of the population mean in Equation 20 with respect to the allele frequency such that
which, similar to the theoretical form in quantitative genetics, also contains both the additive effect of the $j$th SNP and additional terms encoding the interaction effect between the $j$th SNP and all other variants in the data. Once again, under HardyWeinberg equilibrium, the additive variance for the $j$th SNP is found as taking on the following form
where we can explicitly draw connections between the two frameworks by setting ${\beta}_{j}={a}_{j}$ and ${\theta}_{jk}={[aa]}_{jk}$. Note that when there no nonadditive effects (such that $\mathit{\bm{\theta}}=\mathrm{\U0001d7ce}$), the above reduces to ${\sigma}_{A}^{2}={\sum}_{j}2{p}_{j}(1{p}_{j}){\beta}_{j}^{2}$ which resembles the classical form for the additive genetic variance (Lynch and Walsh, 1998).
Full derivation of interaction LD score regression
Request a detailed protocolIn order to derive the interaction LD score (iLDSC) regression framework, recall that our goal is to recover missing heritability from GWAS summary statistics by incorporating an additional score that measures the nonadditive genetic variation that is tagged by genotyped SNPs. To do this, we build upon the LD score regression framework and the LDSC software (BulikSullivan et al., 2015b). Here, we assume nonzero contributions from cisacting pairwise interaction effects in the generative model of complex traits as in Equation 16, and we use the observed least squares estimates from Equation 11 to compute chisquare statistics ${\chi}_{j}^{2}=N{\widehat{\beta}}_{j}^{2}$ for every $j=1,\mathrm{\dots},J$ variant in the data. Taking the expectation of these statistics yields
We can simplify Equation 23 in two steps. First, by combining the prior assumption in Equation 8 and the asymptotic approximation in Equation 16, we can show that marginal expectation (i.e. when not conditioning on the true coefficients) $\mathbb{E}[{\widehat{\beta}}_{j}]=0$ for all variants. Second, by conditioning on the generative model from Equation 7, we can use the law of total variance to simplify $\mathbb{V}[{\widehat{\beta}}_{j}]$ where
since ${\mathbf{\mathbf{x}}}_{j}^{\u22ba}{\mathbf{\mathbf{X}}}_{D}=\mathrm{\U0001d7ce}$. Using the same logic from the original LDSC regression framework (BulikSullivan et al., 2015b), we can use Isserlis’ theorem Isserlis, 1918 to write the above in terms of more familiar quantities based on sample correlations
where ${\stackrel{~}{r}}_{jk}$ is used to denote the sample correlation between additivelycoded genotypes at the $j$th and $k$th variants, and ${\stackrel{~}{v}}_{jm}$ is used to denote the sample correlation between the genotype of the $j$th variant and the $m$th genetic interaction on the phenotype of interest (again see Equation 16). Furthermore, we can use the delta method (only displaying terms up to $\mathcal{O}(1/{N}^{2})$) to show that (in expectation)
Next, we can then approximate the quantities in Equation 24 via the following
where ${\mathrm{\ell}}_{j}$ is the corresponding LD score for the additive effect of the $j$th variant and f_{j} represents the “interaction” LD score between the $j$th SNP and all other variants in the data set (Crawford et al., 2017), respectively. Altogether, this leads to the specification of the univariate framework with the $j$th SNP
where we define $\tau =N{\phi}_{\beta}^{2}/J$ as estimates of the additive genetic signal, the coefficient $\vartheta =N{\phi}_{\theta}^{2}/M$ as an estimate of the proportion of phenotypic variation explained by tagged pairwise interaction effects, and 1 is the intercept meant to model the misestimation due to uncontrolled confounding effects (e.g. cryptic relatedness and population stratification). Similar to the original LDSC formulation, an intercept greater than one means significant bias. Note that the simplification for many of the terms above such as $(1{H}^{2})/N\approx 1/N$ results from our assumption that the number of individuals in our study is large. For example, the sample sizes for each biobankscale study considered in the analyses of this manuscript are at least on the order of $N\ge {10}^{4}$ observations (see Supplementary file 5). Altogether, we can jointly express Equation 27 in multivariate form as
where ${\mathit{\bm{\chi}}}^{2}=({\chi}_{1}^{2},\mathrm{\dots},{\chi}_{J}^{2})$ is a $J$dimensional vector of chisquare summary statistics, and $\mathbf{\ell}=({\mathrm{\ell}}_{1},\mathrm{\dots},{\mathrm{\ell}}_{J})$ and $\mathit{\bm{f}}=({f}_{1},\mathrm{\dots},{f}_{J})$ are $J$dimensional vectors of additive and cisinteraction LD scores, respectively. It is important to note that, while ${\mathit{\bm{\chi}}}^{2}$ must be recomputed for each trait of interest, both vectors $\mathbf{\ell}$ and $\mathit{\bm{f}}$ only need to be constructed once per reference panel or individuallevel genotypes (see next section for efficient computational strategies).
To identify summary statistics that have significant tagged interaction effects, we test the null hypothesis ${H}_{0}:\vartheta =0$. The iLDSC software package implements the same model fitting strategy as LDSC. Here, we use weighted least squares to fit the joint regression in Equation 28 such that
where $\mathbf{\mathbf{\Psi}}$ is a $J\times J$ diagonal weight matrix with nonzero elements set to values inversely proportional to the conditional variance $\mathbb{V}[{\chi}_{j}^{2}{\mathrm{\ell}}_{j},{f}_{j}]={\psi}_{jj}^{1}$ to adjust for both heteroscedasticity and overestimation of the summary statistics for each SNP (BulikSullivan et al., 2015b). Standard errors for each coefficient estimate are derived via a jackknife over blocks of SNPs in the data (Finucane et al., 2015), and we then use those standard errors to derive pvalues with a twosided test (i.e. testing the alternative hypothesis ${H}_{A}:\vartheta \ne 0$). It is worth noting that the blockjackknife approach tends to be conservative and yield larger standard errors for hypothesis testing (Efron, 1982). As an alternative, we could first run iLDSC using the blockjackknife procedure over all traits in a study and then use the average of the standard errors to calculate the statistical significance of coefficient estimates; but we do not explore this strategy here and leave that for future work. The quantitative genetics expression for the additive variance ${\sigma}_{A}^{2}$ in Equation 22 is important because it represents the theoretical upper bound on the proportion of phenotypic variance that can be explained from GWAS summary statistics via iLDSC. Using this relationship, we can write the following (approximate) inequality
For all analyses in this paper, we estimate proportion of phenotypic variance explained by genetic effects using a sum of the coefficients $\widehat{\tau}+\widehat{\vartheta}$ (i.e. the estimated additive component plus the additional genetic variance explained by the tagged pairwise interaction effects).
Efficient computation of cisinteraction LD scores
Request a detailed protocolIn practice, cisinteraction LD scores in iLDSC can be computed efficiently through realizing two key opportunities for optimization. First, given $J$ SNPs, the full matrix of genomewide interaction effects $\mathbf{\mathbf{W}}$ contains on the order of $J(J1)/2$ total pairwise interactions. However, to compute the cisinteraction score for each SNP, we simply can replace the full $\mathbf{\mathbf{W}}$ matrix with a subsetted matrix ${\mathbf{\mathbf{W}}}_{j}$ which includes only interactions involving the $j$th SNP. Analogous to the original LDSC formulation (BulikSullivan et al., 2015b), we consider only interactive SNPs within a ciswindow proximal to the focal $j$th SNP for which we are computing the iLDSC score. In the original LDSC model, this is based on the observation that LD decays outside of a window of 1 centimorgan (cM) (BulikSullivan et al., 2015b); therefore, SNPs outside the 1 cM window centered on the $j$th SNP $j$ will not significantly contribute to its LD score. The second opportunity for optimization comes from the fact that the matrix of interaction effects for any focal SNP, ${\mathbf{\mathbf{W}}}_{j}$, does not need to be explicitly generated. Referencing Equation 24, the iLDSC scores are defined as ${\mathbf{\mathbf{x}}}_{j}^{\u22ba}{\mathbf{\mathbf{W}}}_{j}{\mathbf{\mathbf{W}}}_{j}^{\u22ba}{\mathbf{\mathbf{x}}}_{j}/{N}^{2}$. This can be rewritten as ${\mathbf{\mathbf{x}}}_{j}^{\u22ba}({\mathbf{\mathbf{D}}}_{j}{\mathbf{\mathbf{X}}}^{(j)}){({\mathbf{\mathbf{D}}}_{j}{\mathbf{\mathbf{X}}}^{(j)})}^{\u22ba}{\mathbf{\mathbf{x}}}_{j}$, where ${\mathbf{\mathbf{D}}}_{j}=\text{diag}({\mathbf{\mathbf{x}}}_{j})$ is a diagonal matrix with the $j$th genotype as its nonzero elements (Crawford et al., 2017) and ${\mathbf{\mathbf{X}}}^{(j)}$ denotes the subset SNPs within a ciswindow proximal to the focal $j$th SNP. This means that the iLDSC score for the $j$th SNP can be simply computed as the following
With these simplifications, the computational complexity of generating iLDSC scores reduces to that of computing LD scores — modulo a vectorbyvector Hadamard product which, for each SNP, is constant factor of $N$ (i.e. the number of genotyped individuals).
Coefficient estimates as determined by cisinteraction window size
Request a detailed protocolWhen computing cisinteraction LD scores, the most important decision is choosing the number of interacting SNPs to include in ${\mathbf{\mathbf{X}}}^{(j)}$ (or equivalently ${\mathbf{\mathbf{W}}}_{j}$ for each $j$th focal SNP in the calculation of f_{j} in Equation 31). The iLDSC framework considers different estimating windows to account for our lack of a priori knowledge about the ‘correct’ nonadditive genetic architecture of traits. Theoretically, one could follow previous work Guan and Stephens, 2011; Carbonetto and Stephens, 2012; Zhou et al., 2013; Zhu and Stephens, 2017; Zhu and Stephens, 2018; Demetci et al., 2021 by considering an $L$valued grid of possible SNP interaction window sizes. After fitting a series of iLDSC regressions with cisinteraction LD scores ${\mathit{\bm{f}}}^{(l)}$ generated under the $L$different window sizes, we could compute normalized importance weights using their maximized likelihoods via the following
As a final step in the model fitting procedure, we could then compute averaged estimates of the coefficients $\tau $ and $\vartheta $ by marginalizing (or averaging) over the $L$different grid combinations of estimating windows
This final step can be viewed as an analogy to model averaging where marginal estimates are computed via a weighted average using the importance weights (Hoeting et al., 1999). In the current study, we explore the utility of cisinteraction LD scores generated with different window sizes ± 5, ± 10, ± 25, and ± 50 SNPs around each $j$th focal SNP. In practice, we find that cisinteraction LD scores that are calculated using larger windows lead to the most robust estimates of heritability while also not over representing the total phenotypic variation explained by tagged nonadditive genetic effects (see Figure 3—figure supplement 1). Therefore, unless otherwise stated, we use cisinteraction LD scores calculated with a ± 50 SNP interaction window for all simulations and real data analyses conducted in this work. For a direct comparison between choosing a single window size versus the model averaging strategy described above, see Supplementary files 1 and 2.
Relationship between minor allele frequency and effect size
Request a detailed protocolThe LDSC software computes LD scores using annotations over equally spaced minor allele frequency (MAF) bins. These annotations enable the per trait relationship between the MAF and the effect size of each variant in the genome to vary based on the discrete category (or MAF bin) it is placed into. This additional flexibility is intended to help LDSC be more robust when estimating heritability. The relationship between MAF and effect size is already implicitly encoded in the LDSC formulation since we assume genotypes are normalized. When normalizing by the variance of each SNP (or equivalently its MAF), we make the assumption that rare variants inherently have larger effect sizes. There exists a true functional relationship between MAF and effect size which is likely to be somewhere between the two extremes of (i) normalizing each SNP by its MAF and (ii) allowing the variance per SNP to be dictated by its MAF.
Recent approaches have proposed using a single parameter $\alpha $ to better represent the nonlinear relationship between MAF and variant effect size. The main idea is that this $\alpha $ not only provides the same additional flexibility to LDSC as the MAFbased discrete annotations, but it also empirically yields even more precise heritability estimates (Zabad et al., 2021). Namely, we use
where ${a}_{c}(k)$ is the annotation value for the $c$th categorical bin. The α parameter is unknown in practice and needs to be estimated for any given trait. While standard ranges for α can be used for heritability estimates, we use a restricted maximum likelihood (REML) based method which was recently developed (Schoech et al., 2019).
In the iLDSC software, we use this α construction to handle the relationship between MAF and variant effect size for two specific reasons. First, by constructing the LD scores using α, we more accurately capture the variation in chisquare test statistics due to additive effects (Zabad et al., 2021). Second, we note that there is correlation between MAF and (i) LD scores, (ii) cisinteraction LD scores, and (iii) trait architecture. To that end, if we do not properly condition on MAF, there becomes additional bias, and we may falsely attribute some amount of variation in the chisquare test statistics to LD or the tagged interaction effects. Therefore, in our formulation, we include an α term on the LD scores to condition on this effect. We demonstrate in simulations that this removes the bias introduced by the relationship between MAF and trait architecture, and it mitigates potential inflation of type I error rates in the iLDSC test.
Estimation of allele frequency parameters
Request a detailed protocolIn the main text, we analyzed 25 complex traits in both the UK Biobank and BioBank Japan data sets. In order to account for minor allele frequency (MAF) dependent trait architecture, we calculated $\alpha $ values for each trait that had not been analyzed by previous studies (Schoech et al., 2019). The α estimates for each of the 25 traits analyzed in this study are shown in Supplementary file 4. Intuitively, $\alpha $ parameterizes the weighting of the effects of each individual variant given its frequency in the study cohort and can take on values in the range of [–1,0]. More negative values of $\alpha $ indicate that lower frequency variants contribute more to the observed variation in a trait of interest, whereas values of α closer to zero indicate that common variants contribute a greater amount of variation to observed trait values.
We took α values for 11 traits (again see Supplementary file 4) that had previously been calculated from Schoech et al. For the remaining 14 traits analyzed in this study, we followed the estimation protocol described in the same manuscript. Specifically, using the variants passing the quality control step in our pipeline for 25,000 randomly selected individuals in the UK Biobank cohort, we constructed MAFdependent genetic relatedness matrices for values of $\alpha =\{1,0.95,0.9,\mathrm{\dots},0\}$ using the GRMMAFLD software (Schoech, 2018). We then used the GCTA software (Yang et al., 2011) to obtain heritability and likelihood estimates using REML for each $\alpha $trait pairing. We then fit a traitspecific profile likelihood across the range of α values and estimate the maximum likelihood value of $\alpha $ using a natural cubic spline.
Simulation studies
We used a simulation scheme to generate synthetic quantitative traits and SNPlevel summary statistics under multiple genetic architectures using real genomewide data from individuals of selfidentified European ancestry in the UK Biobank. Here, we consider phenotypes that have some combination of additive effects, cisacting interactions, and a third source of genetic variance stemming from either genebyenvironment (G×E) or genebyancestry (G×Ancestry) effects. For each scenario, we select some set of SNPs to be causal and assume that complex traits are generated via the following general linear model
where $\mathbf{\mathbf{y}}$ is an $N$dimensional vector containing all the phenotypes; $\mathbf{\mathbf{X}}$ is an $N\times J$ matrix of genotypes encoded as 0, 1, or 2 copies of a reference allele; β is a $J$dimensional vector of additive effect sizes for each SNP; $\mathbf{\mathbf{W}}$ is an $N\times M$ matrix which holds all pairwise interactions between the randomly selected subset of the interacting SNPs with corresponding effects θ is an $N\times K$ matrix of either G×E or G×Ancestry interactions with coefficients $\mathit{\bm{\gamma}}$; and $\mathit{\bm{\epsilon}}$ is an $N$dimensional vector of environmental noise. The phenotypic variation is assumed to be $\mathbb{V}[\mathbf{\mathbf{y}}]=1$. All additive and interaction effect sizes for SNPs are randomly drawn from independent standard Gaussian distributions and then rescaled so that they explain a fixed proportion of the phenotypic variance $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]+\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]+\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]={H}^{2}$. Note that we do not assume any specific correlation structure between the effect sizes β, θ, and $\mathit{\bm{\gamma}}$. We then rescale the random error term such that $\mathbb{V}[\mathit{\bm{\epsilon}}]=(1{H}^{2})$. In the main text, we compare the traditional LDSC to its direct extension in iLDSC. For each method, GWAS summary statistics are computed by fitting a singleSNP univariate linear model via least squares where ${\widehat{\beta}}_{j}={({\mathbf{\mathbf{x}}}_{j}^{\u22ba}{\mathbf{\mathbf{x}}}_{j})}^{1}{\mathbf{\mathbf{x}}}_{j}^{\u22ba}\mathbf{\mathbf{y}}$ for every $j=1,\mathrm{\dots},J$ SNP in the data. These effect size estimates are used to derive the chisquare test statistics ${\chi}_{j}^{2}=N{\widehat{\beta}}_{j}^{2}$. We implement both LDSC and iLDSC with the LD matrix $\mathbf{\mathbf{R}}={\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{X}}/N$ and the cisinteraction correlation matrix $\mathbf{\mathbf{V}}={\mathbf{\mathbf{X}}}^{\u22ba}\mathbf{\mathbf{W}}/N$ being computed using a reference panel of 489 individuals from the European superpopulation (EUR) of the 1000 Genomes Project (https://mathgen.stats.ox.ac.uk/impute/data_download_1000G_phase1_integrated.html). The resulting matrices $\mathbf{\mathbf{R}}$ and $\mathbf{\mathbf{V}}$ are used to compute the additive and cisinteraction LD scores, respectively.
Polygenic simulations with cisinteractions
Request a detailed protocolIn our first set of simulations, we consider phenotypes with polygenic architectures that are made up of only additive and cisacting SNPbySNP interactions. Here, we begin by assuming that every SNP in the genome has at least a small additive effect on the traits of interest. Next, when generating synthetic traits, we assume that the additive effects make up $\rho \%$ of the heritability while the pairwise interactions make up the remaining $(1\rho )\%$. Alternatively, the proportion of the heritability explained by additivity is said to be $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=\rho {H}^{2}$, while the proportion detailed by interactions is given as $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=(1\rho ){H}^{2}$. The setting of $\rho =1$ represents the limiting null case for iLDSC where the variation of a trait is driven by solely additive effects. Here, we use the same simulation strategy used in Crawford et al. where we divide the causal cisinteraction variants into two groups. One may view the SNPs in group #1 as being the ‘hubs’ of an interaction map. SNPs in group #2 are selected to be variants within some kilobase (kb) window around each SNP in group #1. Given different parameters for the generative model in Equation 35, we simulate data mirroring a wide range of genetic architectures by toggling the following parameters:
heritability: ${H}^{2}=$ 0.3 and 0.6;
proportion of phenotypic variation that is generated by additive effects: $\rho =$ 0.5, 0.8, and 1;
percentage of SNPs selected to be in group #1: 1% (sparse), 5%, and 10% (polygenic);
genomic window used to assign SNPs to group #2: ± 10 and ± 100 kilobase (kb);
allele frequency parameter: $\alpha =$ −1,–0.5, and 0.
All figures and tables show the mean performances (and standard errors) across 100 simulated replicates.
Polygenic simulations with genebyenvironmental effects
Request a detailed protocolIn our second set of simulations, we continue to consider phenotypes with polygenic architectures that are made up of only additive and cisacting SNPbySNP interactions; however, now we also consider each trait to have contributions stemming from nonzero G×E effects. Here, both the additive and cisinteraction effects are simulated in the same way as previously described where, for the two groups of interacting variants, 10% of SNPs were selected to be in group #1 and we chose ±10 kb windows to assign SNPs to group #2. To create G×E effects, we follow a simulation strategy implemented by Zhu et al. and split our sample population in half to emulate two subsets of individuals coming from different environments. We randomly draw the effect sizes for the first environment from a standard Gaussian distribution which we denote as ${\mathit{\bm{\gamma}}}_{1}$. We then selected an amplification coefficient $w$ and set the effect sizes of the G×E interactions in the second environment to be a scaled version of the first environment effects where ${\mathit{\bm{\gamma}}}_{2}=w{\mathit{\bm{\gamma}}}_{1}$. In this paper, we generate traits with heritability ${H}^{2}=\{0.3,0.6\}$ and amplification coefficients set to $w=[1.1,1.2,\mathrm{\dots},2]$. For the first set of simulations, we hold the proportion of phenotypic variation explained by the different genetic components constant by fixing:
${H}^{2}=0.3$: $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=0.15$; $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=0.075$; and $\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]=0.075$;
${H}^{2}=0.6$: $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=0.3$; $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=0.15$; and $\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]=0.15$;
where $\mathbf{\mathbf{Z}}=[{\mathbf{\mathbf{X}}}_{1},{\mathbf{\mathbf{X}}}_{2}]$ is the set of genotypes split according to environment and $\mathit{\bm{\gamma}}=[{\mathit{\bm{\gamma}}}_{1},{\mathit{\bm{\gamma}}}_{2}]$. To test the sensitivity of the cisinteraction LD scores to other sources of nonadditive variation, we also repeated the same simulations where there were only additive and G×E effects contributing equally to trait architecture:
${H}^{2}=0.3$: $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=0.15$; $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=0$; and $\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]=0.15$;
${H}^{2}=0.6$: $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=0.3$; $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=0$; and $\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]=0.3$.
Again all figures show the mean performances (and standard errors) across 100 simulated replicates.
Polygenic simulations with genebyancestry effects
Request a detailed protocolIn our third set of simulations, we consider phenotypes with polygenic architectures that are made up of additive, cisinteractions, and G×Ancestry effects. Here, we follow Sohail et al. and first run a matrix decomposition on the individuallevel genotype matrix $\mathbf{\mathbf{X}}={\mathbf{\mathbf{U}\mathbf{Q}}}^{\u22ba}$ where $\mathbf{\mathbf{U}}$ is a unitary $N\times K$ score matrix, $\mathbf{\mathbf{Q}}$ is a $K\times J$ loadings matrix, and $K$ represents the number of (predetermined) principal components (PCs). To generate G×Ancestry interactions, we then create the matrix ${\mathbf{\mathbf{Z}}}_{k}={\mathbf{\mathbf{X}\mathbf{q}}}_{k}$ where ${\mathbf{\mathbf{q}}}_{k}$ is a $J$dimensional vector of SNP loadings for the $k$th principal component. In this paper, we generate traits with heritability ${H}^{2}=\{0.3,0.6\}$ and interaction effects taken over $k=1,\mathrm{\dots},10$ principal components. For the first set of simulations, we hold the proportion of phenotypic variation explained by the different genetic components constant by fixing:
${H}^{2}=0.3$: $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=0.15$; $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=0.075$; and $\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]=0.075$;
${H}^{2}=0.6$: $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=0.3$; $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=0.15$; and $\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]=0.15$;
To test the sensitivity of the cisinteraction LD scores to other sources of nonadditive variation, we also repeated the same simulations where there were only additive and G×E effects contributing equally to trait architecture:
${H}^{2}=0.3$: $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=0.15$; $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=0$; and $\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]=0.15$;
${H}^{2}=0.6$: $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]=0.3$; $\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]=0$; and $\mathbb{V}[\mathbf{\mathbf{Z}}\mathit{\bm{\gamma}}]=0.3$.
Note that, for each case, we generate summary statistics in two ways: (i) including the top 10 PCs as covariates in the marginal linear model to correct for population structure and (ii) not correcting for any population structure. Again all figures show the mean performances (and standard errors) across 100 simulated replicates.
Sparse simulation study design with additive effects
Request a detailed protocolIn this set of simulations, we consider phenotypes with sparse architectures (Zhou et al., 2013). Here, traits were simulated with solely additive effects such that $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]={H}^{2}$, but this time only variants with the top or bottom $\{1,5,10,25,50,100\}$ percentile of LD scores were given nonzero coefficients (a similar simulation approach was also previously implemented in both BulikSullivan et al., 2015b and Lee et al., 2018). We once again generate traits with heritability ${H}^{2}=\{0.3,0.6\}$. We also want to note that, in each of these specific analyses, synthetic trait architectures were generated using all UK Biobank genotyped variants that passed initial preprocessing and quality control (see next section). Since not all of these SNPs are HapMap3 SNPs, some variants were omitted from the LDSC and iLDSC regression. Overall, as shown in the main text with results taken over 100 replicates, breaking the assumed relationship between LD scores and chisquared statistics (i.e. that they are generally positively correlated) led to unbounded estimates of heritability in all but the (more polygenic) scenario when 100% of SNPs contributed to phenotypic variation.
Polygenic simulations with unobserved additive effects
Request a detailed protocolIn this next set of simulations, we consider another extension of the polygenic case where a portion of the variants with only additive genetic effects are not observed due ascertainment or other quality control procedures. It was found in Hemani et al., 2014. that an initial set of signals pointing towards evidence of genetic interactions were actually better explained using linear models of unobserved variants in the same haplotype. Here, we test whether the iLDSC framework is prone to overestimate the nonadditive genetic variance when additive effects in the same haplotype are not included in the model. In each simulation, we generated haplotypes that each contain 5000 variants. Next, we select either a single causal variant with only an additive effect or a set of ten causal variants with only additive effects — each having an MAF that is randomly selected between: (i) (0.01, 0.1), (ii) (0.1, 0.2), (iii) (0.2, 0.3), (iv) (0.3, 0.4), and (v) (0.4, 0.5). The corresponding additive effect size for each causal variant across the haplotype is simulated inversely proportional with its MAF. For this analysis, we measure the difference between iLDSC coefficient estimates when every variant is included in the model versus when the haplotype causal variants are omitted for two different trait architectures with broadsense heritability set to ${H}^{2}=$ 0.3 and 0.6. Differences in the component estimates between the observed and unobserved single additive variant models are shown in Figure 3—figure supplement 9A and B. Similar estimates when the larger number of ten additive variants are unobserved in each haplotype are shown in Figure 3—figure supplement 9C and D. If iLDSC was prone to overestimating the nonadditive effects, then the omission of the variants with only significant additive effects would lead to increased estimates of $\tau $ and $\vartheta $. However, across a range of generative broadsense heritabilities and haplotype architectures we observe that estimates of $\tau $ and $\vartheta $ are robust. Intuitively, this is likely due to the fact that these simulations were done under polygenic trait architectures where, as a result, the omission of a few causal variants with small marginal effect sizes has little impact on the ability to estimate genetic variance.
Polygenic simulations with unobserved interaction effects
Request a detailed protocolIn this set of simulations, we extend the polygenic case to a setting where a portion of the variants involved in genetic interactions are unobserved. Similar to the case with unobserved additive effects, the purpose of these simulations is to assess whether the iLDSC framework is prone to false discovery of nonadditive genetic variance when causal interacting SNPs are not included during the estimation of GWAS summary statistics. In each simulation, we generated haplotypes that each contain 5000 variants. Traits were simulated using the generative model in Equation (35) with both additive and interaction effects such that $\mathbb{V}[\mathbf{\mathbf{X}}\mathit{\bm{\beta}}]+\mathbb{V}[\mathbf{\mathbf{W}}\mathit{\bm{\theta}}]={H}^{2}$. Here, every SNP in the genome had at least a small additive effect with a corresponding effect size that was drawn to be inversely proportional to its MAF. Only 1% or 5% of variants within each haplotype had causal nonzero interaction effects. However, when running iLDSC, only a percentage of the interacting SNPs {1%, 5%, 10%, 25% or 50%} were included in the estimation of $\widehat{\vartheta}$. We once again generate traits with heritability ${H}^{2}=\{0.3,0.6\}$ such that the proportion of genetic variance explained by additive effects was equal to $\rho =\{0.5,0.8\}$. As with the other simulation scenarios, all synthetic traits were generated using UK Biobank genotyped variants that passed initial preprocessing and quality control (see next section). Since not all of these SNPs are HapMap3 SNPs, some variants were omitted from the iLDSC regression analyses. Overall, as discussed in the main text with results taken over 100 replicates, iLDSC underestimated values of $\widehat{\vartheta}$ when there were unobserved interacting variants (see Figure 3—figure supplements 10 and 11). As expected, estimates of the additive variance component $\widehat{\tau}$, on the other hand, were not affected.
Polygenic simulations with correlated additive and interaction effects
Request a detailed protocolIn our last set of simulations, we sought out to better understand how the relationship between the additive ($\beta $) and interaction ($\theta $) coefficients in the generative model of complex traits could potentially bias the additive and nonadditive variance component estimates in LDSC and iLDSC. To that end, we performed a set of simulations where we varied the correlation between the set of effects. Specifically, we first drew a set of additive effect sizes for each variant using the MAFdependent procedure described above (i.e. $\alpha =1$). We next selected a subset of the causal variants to be in cisinteractions. Here, we set the interaction effect sizes to covary with the additive effect size vector in two different ways. In the first, we simply drew the additive and interaction effect sizes from a multivariate normal such that their correlation was equal to $r=\{1,0.8,0.6,\mathrm{\dots},0.6,0.8,1\}$ (see Figure 3—figure supplement 12). In the second, we simply amplified the interaction effects to be a linear function $\theta =\beta \times q$ (Figure 3—figure supplement 13A and C) or a squared function $\theta ={\beta}^{2q}$ (Figure 3—figure supplement 13B and D) of the additive effects where $q=\{0.1,0.2,\mathrm{\dots},0.9,1\}$. While testing 100 replicates for each value of $q$, we observed that the mean estimate of genetic variance had a slight upward bias as the correlation between the additive and interaction effect sizes in the generative model increased; however, the distribution of these bias estimates covered zero in the first and third quartiles of all results. We evaluated this behavior for multiple broadsense heritability levels ${H}^{2}$ = 0.3 and 0.6.
Preprocessing for the UK Biobank and BioBank Japan
Request a detailed protocolIn order to apply the iLDSC framework to 25 continuous traits the UK Biobank (Bycroft et al., 2018), we first downloaded genotype data for 488,377 individuals in the UK Biobank using the ukbgene tool (https://biobank.ctsu.ox.ac.uk/crystal/download.cgi) and converted the genotypes using the provided ukbconv tool (https://biobank.ctsu.ox.ac.uk/crystal/refer.cgi?id=149660). Phenotype data for the 25 continuous traits were also downloaded for those same individuals using the ukbgene tool. Individuals identified by the UK Biobank as having high heterozygosity, excessive relatedness, or aneuploidy were removed (1,550 individuals). After separating individuals into selfidentified ancestral cohorts using data field 21000, unrelated individuals were selected by randomly choosing an individual from each pair of related individuals. This resulted in $N=$ 349,469 white British individuals to be included in our analysis. We downloaded imputed SNP data from the UK Biobank for all remaining individuals and removed SNPs with an information score below 0.8. Information scores for each SNP are provided by the UK Biobank (http://biobank.ctsu.ox.ac.uk/crystal/refer.cgi?id=1967).
Quality control for the remaining genotyped and imputed variants was then performed on each cohort separately using the following steps. All structural variants were first removed, leaving only single nucleotide polymorphisms (SNPs) in the genotype data. Next, all AT/CG SNPs were removed to avoid possible confounding due to sequencing errors. Then, SNPs with minor allele frequency less than 1% were removed using the PLINK 2.0 (Chang et al., 2015) command maf 0.01. We then removed all SNPs found to be out of HardyWeinberg equilibrium, using the PLINK hwe 0.000001 flag to remove all SNPs with a Fisher’s exact test pvalue $>{10}^{6}$. Finally, all SNPs with missingness greater than 1% were removed using the PLINK mind 0.01 flag.
We then performed a genomewide association study (GWAS) for each trait in the UK Biobank on the remaining 8,981,412 SNPs. SNPlevel GWAS effect sizes were calculated using PLINK and the glm flag (Chang et al., 2015). Age, sex, and the first 20 principal components were included as covariates for all traits analyzed (Sohail et al., 2019). Principal component analysis was performed using FlashPCA 2.0 (Abraham et al., 2017) on a set of independent markers derived separately for each ancestry cohort using the PLINK command indeppairwise 100 10 0.1. Using the parameters indeppairwise removes all SNPs that have a pairwise correlation above 0.1 within a 100 SNP window, then slides forward in increments of ten SNPs genomewide.
In order to analyze data from BioBank Japan, we downloaded publicly available GWAS summary statistics for the 25 traits listed in Supplementary file 5 from https://pheweb.jp/downloads. Summary statistics used age, sex, and the first ten principal components as confounders in the initial GWAS study. We then used individuals from the East Asian (EAS) superpopulation from the 1000 Genomes Project Phase 3 to calculate paired LDSC and iLDSC scores from a reference panel. We pruned the reference panel using the PLINK command indeppairwise 100 10 0.5 to limit the computational time of calculating scores (Chang et al., 2015). This resulted in reference scores for 1,164,666 SNPs that are included on the iLDSC GitHub repository (https://github.com/lcrawlab/iLDSC). Using summary statistics from BioBank Japan, with scores calculated from the EAS population in the 1000 Genomes, we obtained iLDSC heritability estimates for each of the 25 traits.
Data availability
Source code and tutorials for implementing interactionLD score regression via the iLDSC package are written in Python and are publicly available online on GitHub (copy archived at Crawford and Smith, 2024). Files of LD scores, cisinteraction LD scores, and GWAS summary statistics used for our analyses of the UK Biobank and BioBank Japan can be downloaded from the Harvard Dataverse. All software for the traditional and stratified LD score regression framework with LDSC and sLDSC were fit using the default settings, unless otherwise stated in the main text. Source code for these approaches was downloaded from https://github.com/bulik/ldsc (BulikSullivan et al., 2020). When applying sLDSC, we used 97 functional annotations from Gazal et al., 2017 to estimate heritability. Data from the UK Biobank Resource (Bycroft et al., 2018) was made available under Application Numbers 14649 and 22419. Data can be accessed by direct application to the UK Biobank.

Harvard DataverseReplication Data for: Discovering nonadditive heritability using additive GWAS summary statistics.https://doi.org/10.7910/DVN/W6MA8J
References

BookIntroductory Econometrics: Using Monte Carlo Simulation with Microsoft ExcelCambridge University Press.https://doi.org/10.1017/CBO9780511809231

An atlas of genetic correlations across human diseases and traitsNature Genetics 47:1236–1241.https://doi.org/10.1038/ng.3406

SoftwareInteractionLD score (ILDSC) regression, version swh:1:rev:2d828d50502a341a8148f14cde5825c812a04f90Software Heritage.

Genomic heritability: what is it?PLOS Genetics 11:e1005048.https://doi.org/10.1371/journal.pgen.1005048

Missing heritability and strategies for finding the underlying causes of complex diseaseNature Reviews. Genetics 11:446–450.https://doi.org/10.1038/nrg2809

XV.—The correlation between relatives on the supposition of mendelian inheritanceTransactions of the Royal Society of Edinburgh 52:399–433.https://doi.org/10.1017/S0080456800012163

Bayesian variable selection regression for genomewide association studies and other largescale problemsThe Annals of Applied Statistics 5:1780–1815.https://doi.org/10.1214/11AOAS455

Estimation of nonadditive genetic variance in human complex traits from a large sample of unrelated individualsAmerican Journal of Human Genetics 108:786–798.https://doi.org/10.1016/j.ajhg.2021.02.014

Colocalization of GWAS and eQTL signals detects target genesAmerican Journal of Human Genetics 99:1245–1260.https://doi.org/10.1016/j.ajhg.2016.10.003

Statistical and functional studies identify epistasis of cardiovascular risk genomic variants from genomewide association studiesJournal of the American Heart Association 9:e014146.https://doi.org/10.1161/JAHA.119.014146

Shared heritability of human face and brain shapeNature Genetics 53:830–839.https://doi.org/10.1038/s4158802100827w

Genetic interactions drive heterogeneity in causal variant effect sizes for gene expression and complex traitsAmerican Journal of Human Genetics 109:1286–1297.https://doi.org/10.1016/j.ajhg.2022.05.014

Efficient variance components analysis across millions of genomesNature Communications 11:4020.https://doi.org/10.1038/s41467020175769

Quantifying the contribution of dominance deviation effects to complex trait variation in biobankscale dataAmerican Journal of Human Genetics 108:799–808.https://doi.org/10.1016/j.ajhg.2021.03.018

PLINK: a tool set for wholegenome association and populationbased linkage analysesAmerican Journal of Human Genetics 81:559–575.https://doi.org/10.1086/519795

Contrasting the genetic architecture of 30 complex traits from summary association dataAmerican Journal of Human Genetics 99:139–153.https://doi.org/10.1016/j.ajhg.2016.05.013

Leveraging LD eigenvalue regression to improve the estimation of SNP heritability and confounding inflationAmerican Journal of Human Genetics 109:802–811.https://doi.org/10.1016/j.ajhg.2022.03.013

Allele coding in genomic evaluationGenetics, Selection, Evolution 43:25.https://doi.org/10.1186/129796864325

The influence of higherorder epistasis on biological fitness landscape topographyJournal of Statistical Physics 172:208–225.https://doi.org/10.1007/s1095501819753

Rarevariant association testing for sequencing data with the sequence kernel association testThe American Journal of Human Genetics 89:82–93.https://doi.org/10.1016/j.ajhg.2011.05.029

GCTA: a tool for genomewide complex trait analysisAmerican Journal of Human Genetics 88:76–82.https://doi.org/10.1016/j.ajhg.2010.11.011

Dominance genetic variation contributes little to the missing heritability for human complex traitsAmerican Journal of Human Genetics 96:377–385.https://doi.org/10.1016/j.ajhg.2015.01.001

Bayesian largescale multiple regression with summary statistics from genomewide association studiesThe Annals of Applied Statistics 11:1561–1592.https://doi.org/10.1214/17aoas1046
Article and author information
Author details
Funding
National Institutes of Health (T32 GM128596)
 Samuel Pattillo Smith
 Dana Udwin
National Institutes of Health (RF1 AG073593)
 Samuel Pattillo Smith
National Institutes of Health (R35 GM151108)
 Samuel Pattillo Smith
 Arbel Harpak
National Science Foundation (DMS1439786)
 Gregory Darnell
Alfred P. Sloan Foundation
 Lorin Crawford
David and Lucile Packard Foundation
 Lorin Crawford
National Institutes of Health (R01 GM118652)
 Sohini Ramachandran
National Institutes of Health (R35 GM139628)
 Sohini Ramachandran
National Science Foundation (DBI1452622)
 Sohini Ramachandran
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank Jeffrey P Spence, Roshni Patel, Matthew Aguirre, Mineto Ota, and our anonymous referees for insightful comments on an earlier version of this manuscript as well as the Harpak, Ramachandran, and Crawford Labs for helpful discussions. This research was conducted in part using computational resources and services at the Center for Computation and Visualization at Brown University. This research was also conducted using the UK Biobank Resource under Application Numbers 14649 (LC) and 22419 (SR). SP Smith and D Udwin were trainees supported under the Brown University Predoctoral Training Program in Biological Data Science (NIH T32 GM128596). SP Smith was also supported by NIH RF1AG073593. SP Smith and A Harpak were also supported by NIH R35 GM151108 to A Harpak. G Darnell was supported by NSF Grant No. DMS1439786 while in residence at the Institute for Computational and Experimental Research in Mathematics (ICERM) in Providence, RI. This research was supported in part by an Alfred P Sloan Research Fellowship and a David & Lucile Packard Fellowship for Science and Engineering awarded to L Crawford. This research was also partly supported by US National Institutes of Health (NIH) grant R01 GM118652, NIH grant R35 GM139628, and National Science Foundation (NSF) CAREER award DBI1452622 to S Ramachandran. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of any of the funders.
Copyright
© 2024, Pattillo Smith, Darnell et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 840
 views

 68
 downloads

 0
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Genetics and Genomics
Endocrine disrupting chemicals (EDCs) such as bisphenol S (BPS) are xenobiotic compounds that can disrupt endocrine signaling due to steric similarities to endogenous hormones. EDCs have been shown to induce disruptions in normal epigenetic programming (epimutations) and differentially expressed genes (DEGs) that predispose disease states. Most interestingly, the prevalence of epimutations following exposure to many EDCs persists over multiple generations. Many studies have described direct and prolonged effects of EDC exposure in animal models, but many questions remain about molecular mechanisms by which EDCinduced epimutations are introduced or subsequently propagated, whether there are cell typespecific susceptibilities to the same EDC, and whether this correlates with differential expression of relevant hormone receptors. We exposed cultured pluripotent (iPS), somatic (Sertoli and granulosa), and primordial germ celllike (PGCLC) cells to BPS and found that differential incidences of BPSinduced epimutations and DEGs correlated with differential expression of relevant hormone receptors inducing epimutations near relevant hormone response elements in somatic and pluripotent, but not germ cell types. Most interestingly, we found that when iPS cells were exposed to BPS and then induced to differentiate into PGCLCs, the prevalence of epimutations and DEGs was largely retained, however, >90% of the specific epimutations and DEGs were replaced by novel epimutations and DEGs. These results suggest a unique mechanism by which an EDCinduced epimutated state may be propagated transgenerationally.

 Genetics and Genomics
Microbial secondary metabolites are a rich source for pharmaceutical discoveries and play crucial ecological functions. While tools exist to identify secondary metabolite clusters in genomes, precise sequencetofunction mapping remains challenging because neither function nor substrate specificity of biosynthesis enzymes can accurately be predicted. Here, we developed a knowledgeguided bioinformatic pipeline to solve these issues. We analyzed 1928 genomes of Pseudomonas bacteria and focused on ironscavenging pyoverdines as model metabolites. Our pipeline predicted 188 chemically different pyoverdines with nearly 100% structural accuracy and the presence of 94 distinct receptor groups required for the uptake of ironloaded pyoverdines. Our pipeline unveils an enormous yet overlooked diversity of siderophores (151 new structures) and receptors (91 new groups). Our approach, combining feature sequence with phylogenetic approaches, is extendable to other metabolites and microbial genera, and thus emerges as powerful tool to reconstruct bacterial secondary metabolism pathways based on sequence data.