Discovering non-additive heritability using additive GWAS summary statistics
Abstract
LD score regression (LDSC) is a method to estimate narrow-sense heritability from genome-wide association study (GWAS) summary statistics alone, making it a fast and popular approach. In this work, we present interaction-LD score (i-LDSC) 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 re-analyzing 25 well-studied 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 cis-interaction 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, i-LDSC detects additional variation contributed by genetic interactions. The i-LDSC software and its application to these biobanks represent a step towards resolving further genetic contributions of sources of non-additive 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 genome-wide association summary statistics. The authors present a convincing method, i-LDSC, 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 (Bulik-Sullivan et al., 2015b, Bulik-Sullivan 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 biobank-scale data sets of unrelated individuals (Bulik-Sullivan et al., 2015b; Shi et al., 2016; Hou et al., 2019; Pazokitoroudi et al., 2020). Computational and privacy considerations with genome-wide 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 (Bulik-Sullivan 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 (Bulik-Sullivan et al., 2015b), which corrects for inflation in GWAS summary statistics by modeling the relationship between the variance of SNP-level 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 chi-square 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 Bulik-Sullivan et al. showed that violations of these assumptions can be tolerated to a point, but begin to affect the estimation of narrow-sense 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 (Bulik-Sullivan 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 off-diagonal elements. This tradeoff helps LDSC to scale genome-wide, 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 gene-by-gene (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, non-additive 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 SNP-level GWAS summary statistics can provide evidence of non-additive genetic effects contributing to trait architecture if there is a nonzero correlation between individual-level genotypes and their statistical interactions. We present the ‘interaction-LD score’ regression model or i-LDSC: 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 i-LDSC in simulation studies, we focus on synthetic trait architectures that have been generated with contributions stemming from second-order and cis-acting statistical SNP-by-SNP interaction effects; however, note that the general concept underlying i-LDSC can easily be extended to other sources of non-additive genetic effects (e.g. gene-by-environment interactions). The main difference between i-LDSC and LDSC is that the i-LDSC model includes an additional set of ‘cis-interaction’ 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 i-LDSC recovers substantial non-additive heritability that is not captured by LDSC when genetic interactions are indeed present in the generative model for a given complex trait. More importantly, i-LDSC 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 non-additive 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 i-LDSC, with our development of a new cis-interaction score, represents a significant step towards resolving the true contribution of genetic interactions.
Results
Overview of the interaction-LD score regression model
Interaction-LD score regression (i-LDSC) 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 i-LDSC method and its corresponding software, as well as detail how its underlying model differs from that of LDSC (Bulik-Sullivan et al., 2015b). We will assume that we are analyzing a GWAS dats set where is an matrix of genotypes with denoting the number of SNPs (each of which is encoded as {0, 1, 2} copies of a reference allele at each locus ) and is an -dimensional vector of measurements of a quantitative trait. The i-LDSC framework only requires summary statistics of individual-level data: namely, marginal effect size estimates for each SNP and a sample LD matrix (which can be provided via reference panel data).
We begin by considering the following generative linear model for complex traits
where is an intercept term; is a -dimensional vector containing the true additive effect sizes for an additional copy of the reference allele at each locus on ; is an matrix of (pairwise) cis-acting SNP-by-SNP statistical interactions between some subset of causal SNPs, where columns of this matrix are assumed to be the Hadamard (element-wise) product between genotypic vectors of the form for the -th and -th variants; is an -dimensional vector containing the interaction effect sizes; is a normally distributed error term with mean zero and variance scaled according to the proportion of phenotypic variation not explained by genetic effects (Bulik-Sullivan et al., 2015b), which we will refer to as the broad-sense heritability of the trait denoted by ; and denotes an identity matrix. For convenience, we will assume that the genotype matrix (column-wise) and the trait of interest have been mean-centered and standardized (Strandén and Christensen, 2011; de Los Campos et al., 2013; Zhou et al., 2013). Lastly, we will let the intercept term b0 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 is the proportion of phenotypic variation contributed by additive SNP effects under the generative model, while 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 i-LDSC accurately recovers non-additive 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 non-additive components in Equation 1 are orthogonal to each other. In other words, . 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 and the matrix of genetic interactions 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 (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 non-additive 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 and are not completely orthogonal (i.e. such that ) and the following relationship between the moment matrix , the observed marginal GWAS summary statistics , and the true coefficient values from the generative model in Equation 1 holds in expectation (see Materials and methods)
where is a sample estimate of the LD matrix, and represents a sample estimate of the correlation between the individual-level genotypes and the span of genetic interactions between causal SNPs in . Intuitively, the term 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) non-additive genetic effects do not contribute to the overall architecture of a trait (i.e. such that ) or (ii) the genotype and interaction matrices and 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 i-LDSC is to increase estimates of genetic variance by accounting for sources of non-additive 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 (Bulik-Sullivan et al., 2015b). Here, according to Equation 3, we note that where is a scale variance term due to uncontrolled confounding effects (Guan and Stephens, 2011; Song et al., 2022). Next, we condition on and take the expectation of chi-square statistics to yield
We define as the LD score for the additive effect of the -th variant (Bulik-Sullivan et al., 2015b), and represents the ‘cis-interaction’ LD score which encodes the pairwise interaction between the -th variant and all other variants within a genomic window that is a pre-specified 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 (Bulik-Sullivan et al., 2015b; Song et al., 2022), we get the following simplified regression model
where is a -dimensional vector of chi-square summary statistics, and and are -dimensional vectors of additive and cis-interaction LD scores, respectively. Furthermore, we define the variance components and as the additive and non-additive 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 cis-interaction LD scores by considering only a subset of interactions between each -th focal SNP and SNPs within a cis-proximal window around the -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 -th SNP will not significantly contribute to its LD scores. Note that the width of this window can be relaxed in the i-LDSC software when appropriate. We fit the i-LDSC model using weighted least squares to estimate regression parameters and derive p-values for identifying traits that have significant statistical evidence of tagged cis-interaction effects by testing the null hypothesis . Importantly, under the null model of a trait being generated by only additive effects, the i-LDSC 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 () and interaction ) LD scores are lowly correlated. This is important because it indicates that the presence of cis-interaction LD scores in the model specified in Equation 5 has little-to-no influence over the estimate for the additive coefficient . Instead, the inclusion of creates a multivariate model that can identify the proportion of variance explained by both additive and non-additive effects in summary statistics. In other words, we can interpret as an estimate of the phenotypic variation explained by tagged cis-acting interaction effects. The concept of additive genetic effects partially explaining non-additive variation has also described in various studies from quantitative genetics (Hill et al., 2008; Hivert et al., 2021; Mäki-Tanila and Hill, 2014). Under Hardy-Weinberg equilibrium, it can be shown that the additive variance explained by SNPs takes on the following form (Materials and methods) (Falconer and Mackay, 1983)
The expression for the additive variance 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 i-LDSC framework. As a result, we use the sum of coefficient estimates to construct i-LDSC heritability estimates. A full derivation of the cis-interaction regression framework and details about its corresponding implementation in our software i-LDSC can be found in Materials and Methods.
Detection of tagged pairwise interaction effects using i-LDSC in simulations
We illustrate the power of i-LDSC across different genetic trait architectures via extensive simulation studies (Materials and methods). We generate synthetic phenotypes using real genome-wide genotype data from individuals of self-identified 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 cis-interaction 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: H2 = 0.3 and 0.6;
proportion of phenotypic variation that is generated by additive effects: 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 i-LDSC robustly detects significant tagged non-additive genetic variance, regardless of the total number of causal interactions genome-wide. Instead, the power of i-LDSC 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 cis-interaction LD scores (see Materials and methods). i-LDSC shows a similar performance for detecting tagged cis-interaction 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, i-LDSC does not falsely identify putative non-additive genetic effects in GWAS summary statistics when the synthetic phenotype was generated by only additive effects (). Figure 2 illustrates the performance of i-LDSC under the null hypothesis , with the type I error rates for different estimation window sizes of the cis-interaction LD scores highlighted in panel A. Here, we also show that, when no genetic interaction effects are present, i-LDSC unbiasedly estimates the cis-interaction coefficient in the regression model to be (Figure 2B), robustly estimates the heritability (Figure 2C), and provides well-calibrated p-values when assessed over many traits (Figure 2D). This behavior is consistent across different MAF-dependent effect size distributions, and p-value calibration is not sensitive to misspecification of the estimation windows used to generate the cis-interaction LD scores (Figure 2—figure supplements 1–2).
One of the innovations that i-LDSC offers over the traditional LDSC framework is increased heritability estimates after the identification of non-additive 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 cis-interaction 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 i-LDSC can yield upwardly biased heritability estimates when the cis-interaction 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, i-LDSC accurately partitions genetic effects into additive versus cis-interacting components, which in turn generally leads the ability of i-LDSC to capture more genetic variance. The mean absolute error between the true generative heritability and heritability estimates produced by i-LDSC 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 i-LDSC 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 non-additive genetic effects. Specifically, we compare heritability estimates from LDSC and i-LDSC in the presence of additive effects, cis-acting interactions, and a third source of genetic variance stemming from either gene-by-environment (G×E) or or gene-by-ancestry (G×Ancestry) effects. Details on how these components were generated can be found in Materials and Methods. In general, i-LDSC underestimates overall heritability when additive effects and cis-acting 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 i-LDSC can yield unbounded heritability estimates greater than 1 (Figure 3—figure supplement 4). Also interestingly, when we omit cis-interactions 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), i-LDSC will still estimate a nonzero genetic variance component with the cis-interaction LD scores (Figure 3—figure supplements 5–7). Collectively, these results empirically show the important point that cis-interaction 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) 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 cis-interactions 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 i-LDSC help recover the relationship between the expectation of test statistics and LD.
As a further demonstration of how i-LDSC performs when assumptions of the original LD score model are violated, we also generated synthetic phenotypes with sparse architectures using the spike-and-slab model (Zhou et al., 2013). Here, traits were simulated with solely additive effects, but this time only variants with the top or bottom percentile of LD scores were given nonzero effects (see Materials and methods). Breaking the relationship assumed under the LDSC framework between LD scores and chi-squared 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 i-LDSC estimates of non-additive 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 i-LDSC (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 i-LDSC 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 i-LDSC; 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 i-LDSC framework protects against the false discovery of non-additive genetic effects and underestimates the variance component 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 in i-LDSC 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 i-LDSC 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 i-LDSC 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 i-LDSC 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 pre-computed summary statistics for BioBank Japan were downloaded directly from the consortium website (https://pheweb.jp/downloads). We release the cis-acting SNP-by-SNP interaction LD scores used in our analyses on the i-LDSC 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 cis-interactions (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 non-additive 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 non-additive 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 i-LDSC does (Table 1). In BioBank Japan, 23 of the 25 traits analyzed had a significant nonlinear component detected by i-LDSC — with HDL and triglyceride levels being the only exceptions.
For each of the 25 traits that we analyzed, we found that the i-LDSC heritability estimates are significantly correlated with corresponding estimates from LDSC in both the UK Biobank (, ) and BioBank Japan (, ) 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 (, ) and i-LDSC (, ) analyses as shown in Figure 4B. After comparing the i-LDSC heritability estimates to LDSC, we then assessed whether there was significant difference in the amount of phenotypic variation explained by the non-additive genetic effect component in the GWAS summary statistics derived from the the UK Biobank and BioBank Japan (i.e. comparing the estimates of ; 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 (, ). Notably, the trait with the most significant evidence of tagged cis-interaction effects in GWAS summary statistics is height which is known to have a highly polygenic architecture.
The intercepts estimated by LDSC and i-LDSC 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 i-LDSC intercept estimates tend to be correlated with, but are generally different than, those computed with LDSC — empirically indicating that non-additive 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 (, ) and BioBank Japan (, ).
Lastly, we performed an additional analysis in the UK Biobank where the cis-interaction scores are included as an annotation alongside 97 other functional categories in the stratified-LD score regression framework and its software s-LDSC (Gazal et al., 2017; Materials and methods). Here, s-LDSC 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 () to 0.055 for HDL (; Figure 4C and D). Furthermore, in this analysis, the estimates of the non-additive 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 i-LDSC framework to identify sources of ‘missing’ phenotypic variance explained in heritability estimation. Importantly, moving forward, we suggest using the cis-interaction scores with additional annotations whenever they are available as it provides more conservative estimates of the role of non-additive effects on trait architecture.
Discussion
In this paper, we present i-LDSC, 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 non-additive genetic variation that is tagged by genotyped SNPs. Here, we demonstrate how i-LDSC builds upon the original LDSC model through the development of new ‘cis-interaction’ LD scores which help to investigate signals of cis-acting SNP-by-SNP interactions (Figure 1 and Figure 1—figure supplements 1–5). Through extensive simulations, we show that i-LDSC is well-calibrated 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 i-LDSC provides greater heritability estimates over LDSC when traits are indeed generated with cis-acting SNP-by-SNP interaction effects (Figure 3 and Figure 3—figure supplement 1, and Supplementary files 1 and 2), and we tested the robustness of i-LDSC 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 cis-interaction effects in the UK Biobank and BioBank Japan (Figure 4 and Figure 4—figure supplement 1, Tables 1 and 2, and Supplementary files 3-5). We have made i-LDSC 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 pre-computed cis-interaction 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 i-LDSC framework offers many directions for future development and applications. First, an area of future work would be to explore how the relationship between cis-interaction LD scores and interaction effect sizes from the generative model of complex traits might bias heritability estimates provided by i-LDSC (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 cis-interaction LD scores that we propose are not always enough to recover explainable non-additive genetic effects for all types of trait architectures. While we focus on pairwise cis-acting SNP-by-SNP statistical interactions in this work, the theoretical concepts underlying i-LDSC 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 non-additive 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 i-LDSC 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 genome-wide 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 denotes an -dimensional vector of phenotypic states for a quantitative trait of interest measured in individuals; is an intercept term; is an matrix of genotypes, with denoting the number of single nucleotide polymorphism (SNPs) encoded as copies of a reference allele at each locus; is a -dimensional vector containing the true additive effect sizes for an additional copy of the reference allele at each locus on ; is an matrix that represents the dominance for each genotype encoded as with corresponding effect sizes ; is an matrix of genetic interactions; is an -dimensional vector containing the interaction effect sizes; is a normally distributed error term with mean zero and variance scaled according to the proportion of phenotypic variation not explained by the broad-sense heritability of the trait, denoted by ; and denotes an identity matrix. Note that the encoding for dominance in was chosen because it imposes orthogonality with the genotype encoding in (Purcell et al., 2007; Vitezica et al., 2017; Palmer et al., 2023).
For convenience, we will assume that the genotype matrix (column-wise), the dominance matrix (also column-wise), 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 could encode any source of non-additive genetic interactions (e.g. gene-by-environmental effects) in theory, we limit our focus in this study to trait architectures that have been generated with contributions stemming from cis-acting statistical SNP-by-SNP (or pairwise) interactions. To that end, we assume that the columns of are the Hadamard (element-wise) product between genotypic vectors of the form for the -th and -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 non-independence due to other environmental factors.
As a final set of assumptions, we will let the intercept term 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 and . The broad-sense heritability of the trait is defined as . Under the generative model in Equation 7, we then say that is the proportion of phenotypic variation contributed by additive SNP effects, 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 . 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 i-LDSC accurately recovers non-additive genetic variance in Equation 7 under a broad range of conditions.
Orthogonality between additive and non-additive genetic effects
Request a detailed protocolAssuming that the effect sizes in Equation 8 follow independent and zero mean Gaussian distributions leads to orthogonality between the additive and non-additive components in Equation 7. Since the genotypes and the dominance values are fixed orthogonal matrices, it is straightforward to show that (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 and denoting the -th and -th column of the individual-level genotype matrix and the interaction matrix , 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 and 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 (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 for some constant . However, despite their linear independence, and are themselves not orthogonal and still have a nonzero correlation. This implies that the inner product between genotypes and their interactions is nonzero . To see this, we focus on a focal SNP and consider three different types of interactions:
Scenario I: Interaction between a focal SNP with itself ().
Scenario II: Interaction between a focal SNP with a different SNP ().
Scenario III: Interaction between a focal SNP with a pair of different SNPs ().
The following derivations rely on the fact that: (1) we assume that genotypes have been mean-centered and scaled to have unit variance, and (2) under Hardy-Weinberg equilibrium, SNPs marginally follow a binomial distribution where 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 . With mean-centered SNPs, this is proportional to which is the skewness of the binomial distribution where, again, MAF and 1-MAF of the -th SNP.
Scenario II
Request a detailed protocolAssume that we have two SNPs, and where and represent their respective minor allele frequencies. We want to compute the correlation between and the interaction where . Again, with the mean-centered assumption, the covariance is proportional to the expectation . Here, with SNPs taking on values , the joint distribution between and can be written out as the following Kang and Jung, 2001:
where and is the Pearson correlation or linkage disequilibrium (LD) between the -th and -th SNPs.
Scenario III
Request a detailed protocolThe covariance between a focal SNP and an interaction with a pair of different SNPs will be nonzero if the -th SNP is correlated with either variant (i.e., or ).
Traditional estimation of additive GWAS summary statistics
Request a detailed protocolAs previously mentioned, the key to this work is that SNP-level GWAS summary statistics can also tag non-additive genetic effects when there is a nonzero correlation between individual-level genotypes and their interactions (as defined in Equation 7). Throughout the rest of this section, we will use to denote the LD or pairwise correlation matrix between SNPs. We will then let 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 is again defined as the Pearson correlation coefficient between the -th and -th SNPs, respectively.
In traditional GWAS studies, summary statistics of the true additive effects 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 and the linear effect size estimates :
The second key point combines Equations 10 and 12 to describe the asymptotic relationship between the observed marginal GWAS summary statistics and the joint coefficient values where (in expectation)
After some algebra, the above mirrors a high-dimensional regression model (in expectation) where 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 high-dimensional 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. ), the marginal GWAS summary statistics derived using least squares in Equation 11 will also explain non-additive 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 but does not account for contributions from the non-additive 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 . This simplifies the above to be the following
where the matrix (which we showed to be nonzero) can be interpreted as the sample correlation between individual-level genotypes and the cis-interactions 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 and the true coefficient values
which results from our initial assumption that the residuals are normally distributed with mean zero in Equation 7. Here, we define to represent a sample estimate of the correlation between the individual-level genotypes and the non-additive genetic interaction matrix such that . Similar to the LD matrix , the correlation matrix is also assumed to be computed from reference panel data. Intuitively, when there is additional phenotypic variation contributed by pairwise interactions that can be explained by GWAS effect size estimates. Moreover, when , 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 non-additive variation has also described in classical quantitative genetics (Hill et al., 2008; Hivert et al., 2021; Mäki-Tanila and Hill, 2014). Consider an individual genotyped at loci each with major and minor alleles A and B, respectively. Let be the allele frequency of A at the -th locus, denote the additive effect, and be the additive-by-additive (pairwise) interaction effect between loci and , and represent a third order interaction between loci , , and . 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 additive-by-additive 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 Hardy-Weinberg 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äki-Tanila and Hill, 2014; Hivert et al., 2021). For the -th SNP, the average effect (including terms up to second-order interaction) is given by the following
which notably contains both the additive effect and a summation of additive-by-additive interactions between pairs of loci. The additive genetic variance for the -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 -th locus (again assuming that SNPs marginally follow a binomial distribution). The total additive variance is then obtained by summing over the loci such that (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 individuals, . Ignoring the assumed small contributions from dominance effects, the population mean for a quantitative trait can be written as the following
To find the average effect for the -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 -th SNP and additional terms encoding the interaction effect between the -th SNP and all other variants in the data. Once again, under Hardy-Weinberg equilibrium, the additive variance for the -th SNP is found as taking on the following form
where we can explicitly draw connections between the two frameworks by setting and . Note that when there no non-additive effects (such that ), the above reduces to 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 (i-LDSC) regression framework, recall that our goal is to recover missing heritability from GWAS summary statistics by incorporating an additional score that measures the non-additive genetic variation that is tagged by genotyped SNPs. To do this, we build upon the LD score regression framework and the LDSC software (Bulik-Sullivan et al., 2015b). Here, we assume nonzero contributions from cis-acting 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 chi-square statistics for every 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) for all variants. Second, by conditioning on the generative model from Equation 7, we can use the law of total variance to simplify where
since . Using the same logic from the original LDSC regression framework (Bulik-Sullivan 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 is used to denote the sample correlation between additively-coded genotypes at the -th and -th variants, and is used to denote the sample correlation between the genotype of the -th variant and the -th genetic interaction on the phenotype of interest (again see Equation 16). Furthermore, we can use the delta method (only displaying terms up to ) to show that (in expectation)
Next, we can then approximate the quantities in Equation 24 via the following
where is the corresponding LD score for the additive effect of the -th variant and fj represents the “interaction” LD score between the -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 -th SNP
where we define as estimates of the additive genetic signal, the coefficient 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 results from our assumption that the number of individuals in our study is large. For example, the sample sizes for each biobank-scale study considered in the analyses of this manuscript are at least on the order of observations (see Supplementary file 5). Altogether, we can jointly express Equation 27 in multivariate form as
where is a -dimensional vector of chi-square summary statistics, and and are -dimensional vectors of additive and cis-interaction LD scores, respectively. It is important to note that, while must be recomputed for each trait of interest, both vectors and only need to be constructed once per reference panel or individual-level genotypes (see next section for efficient computational strategies).
To identify summary statistics that have significant tagged interaction effects, we test the null hypothesis . The i-LDSC 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 is a diagonal weight matrix with nonzero elements set to values inversely proportional to the conditional variance to adjust for both heteroscedasticity and over-estimation of the summary statistics for each SNP (Bulik-Sullivan 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 p-values with a two-sided test (i.e. testing the alternative hypothesis ). It is worth noting that the block-jackknife approach tends to be conservative and yield larger standard errors for hypothesis testing (Efron, 1982). As an alternative, we could first run i-LDSC using the block-jackknife 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 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 i-LDSC. 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 (i.e. the estimated additive component plus the additional genetic variance explained by the tagged pairwise interaction effects).
Efficient computation of cis-interaction LD scores
Request a detailed protocolIn practice, cis-interaction LD scores in i-LDSC can be computed efficiently through realizing two key opportunities for optimization. First, given SNPs, the full matrix of genome-wide interaction effects contains on the order of total pairwise interactions. However, to compute the cis-interaction score for each SNP, we simply can replace the full matrix with a subsetted matrix which includes only interactions involving the -th SNP. Analogous to the original LDSC formulation (Bulik-Sullivan et al., 2015b), we consider only interactive SNPs within a cis-window proximal to the focal -th SNP for which we are computing the i-LDSC score. In the original LDSC model, this is based on the observation that LD decays outside of a window of 1 centimorgan (cM) (Bulik-Sullivan et al., 2015b); therefore, SNPs outside the 1 cM window centered on the -th SNP 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, , does not need to be explicitly generated. Referencing Equation 24, the i-LDSC scores are defined as . This can be re-written as , where is a diagonal matrix with the -th genotype as its nonzero elements (Crawford et al., 2017) and denotes the subset SNPs within a cis-window proximal to the focal -th SNP. This means that the i-LDSC score for the -th SNP can be simply computed as the following
With these simplifications, the computational complexity of generating i-LDSC scores reduces to that of computing LD scores — modulo a vector-by-vector Hadamard product which, for each SNP, is constant factor of (i.e. the number of genotyped individuals).
Coefficient estimates as determined by cis-interaction window size
Request a detailed protocolWhen computing cis-interaction LD scores, the most important decision is choosing the number of interacting SNPs to include in (or equivalently for each -th focal SNP in the calculation of fj in Equation 31). The i-LDSC framework considers different estimating windows to account for our lack of a priori knowledge about the ‘correct’ non-additive 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 -valued grid of possible SNP interaction window sizes. After fitting a series of i-LDSC regressions with cis-interaction LD scores generated under the -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 and by marginalizing (or averaging) over the -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 cis-interaction LD scores generated with different window sizes ± 5, ± 10, ± 25, and ± 50 SNPs around each -th focal SNP. In practice, we find that cis-interaction 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 non-additive genetic effects (see Figure 3—figure supplement 1). Therefore, unless otherwise stated, we use cis-interaction 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 to better represent the nonlinear relationship between MAF and variant effect size. The main idea is that this not only provides the same additional flexibility to LDSC as the MAF-based discrete annotations, but it also empirically yields even more precise heritability estimates (Zabad et al., 2021). Namely, we use
where is the annotation value for the -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 i-LDSC 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 chi-square test statistics due to additive effects (Zabad et al., 2021). Second, we note that there is correlation between MAF and (i) LD scores, (ii) cis-interaction 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 chi-square 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 i-LDSC 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 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, 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 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 MAF-dependent genetic relatedness matrices for values of using the GRM-MAF-LD software (Schoech, 2018). We then used the GCTA software (Yang et al., 2011) to obtain heritability and likelihood estimates using REML for each -trait pairing. We then fit a trait-specific profile likelihood across the range of α values and estimate the maximum likelihood value of using a natural cubic spline.
Simulation studies
We used a simulation scheme to generate synthetic quantitative traits and SNP-level summary statistics under multiple genetic architectures using real genome-wide data from individuals of self-identified European ancestry in the UK Biobank. Here, we consider phenotypes that have some combination of additive effects, cis-acting interactions, and a third source of genetic variance stemming from either gene-by-environment (G×E) or gene-by-ancestry (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 is an -dimensional vector containing all the phenotypes; is an matrix of genotypes encoded as 0, 1, or 2 copies of a reference allele; β is a -dimensional vector of additive effect sizes for each SNP; is an matrix which holds all pairwise interactions between the randomly selected subset of the interacting SNPs with corresponding effects θ is an matrix of either G×E or G×Ancestry interactions with coefficients ; and is an -dimensional vector of environmental noise. The phenotypic variation is assumed to be . 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 . Note that we do not assume any specific correlation structure between the effect sizes β, θ, and . We then rescale the random error term such that . In the main text, we compare the traditional LDSC to its direct extension in i-LDSC. For each method, GWAS summary statistics are computed by fitting a single-SNP univariate linear model via least squares where for every SNP in the data. These effect size estimates are used to derive the chi-square test statistics . We implement both LDSC and i-LDSC with the LD matrix and the cis-interaction correlation matrix 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 and are used to compute the additive and cis-interaction LD scores, respectively.
Polygenic simulations with cis-interactions
Request a detailed protocolIn our first set of simulations, we consider phenotypes with polygenic architectures that are made up of only additive and cis-acting SNP-by-SNP 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 of the heritability while the pairwise interactions make up the remaining . Alternatively, the proportion of the heritability explained by additivity is said to be , while the proportion detailed by interactions is given as . The setting of represents the limiting null case for i-LDSC 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 cis-interaction 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: 0.3 and 0.6;
proportion of phenotypic variation that is generated by additive effects: 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: −1,–0.5, and 0.
All figures and tables show the mean performances (and standard errors) across 100 simulated replicates.
Polygenic simulations with gene-by-environmental 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 cis-acting SNP-by-SNP interactions; however, now we also consider each trait to have contributions stemming from nonzero G×E effects. Here, both the additive and cis-interaction 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 . We then selected an amplification coefficient 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 . In this paper, we generate traits with heritability and amplification coefficients set to . For the first set of simulations, we hold the proportion of phenotypic variation explained by the different genetic components constant by fixing:
: ; ; and ;
: ; ; and ;
where is the set of genotypes split according to environment and . To test the sensitivity of the cis-interaction LD scores to other sources of non-additive variation, we also repeated the same simulations where there were only additive and G×E effects contributing equally to trait architecture:
: ; ; and ;
: ; ; and .
Again all figures show the mean performances (and standard errors) across 100 simulated replicates.
Polygenic simulations with gene-by-ancestry effects
Request a detailed protocolIn our third set of simulations, we consider phenotypes with polygenic architectures that are made up of additive, cis-interactions, and G×Ancestry effects. Here, we follow Sohail et al. and first run a matrix decomposition on the individual-level genotype matrix where is a unitary score matrix, is a loadings matrix, and represents the number of (predetermined) principal components (PCs). To generate G×Ancestry interactions, we then create the matrix where is a -dimensional vector of SNP loadings for the -th principal component. In this paper, we generate traits with heritability and interaction effects taken over principal components. For the first set of simulations, we hold the proportion of phenotypic variation explained by the different genetic components constant by fixing:
: ; ; and ;
: ; ; and ;
To test the sensitivity of the cis-interaction LD scores to other sources of non-additive variation, we also repeated the same simulations where there were only additive and G×E effects contributing equally to trait architecture:
: ; ; and ;
: ; ; and .
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 , but this time only variants with the top or bottom percentile of LD scores were given nonzero coefficients (a similar simulation approach was also previously implemented in both Bulik-Sullivan et al., 2015b and Lee et al., 2018). We once again generate traits with heritability . 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 i-LDSC regression. Overall, as shown in the main text with results taken over 100 replicates, breaking the assumed relationship between LD scores and chi-squared 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 i-LDSC framework is prone to overestimate the non-additive 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 i-LDSC coefficient estimates when every variant is included in the model versus when the haplotype causal variants are omitted for two different trait architectures with broad-sense heritability set to 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 i-LDSC was prone to overestimating the non-additive effects, then the omission of the variants with only significant additive effects would lead to increased estimates of and . However, across a range of generative broad-sense heritabilities and haplotype architectures we observe that estimates of and 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 i-LDSC framework is prone to false discovery of non-additive 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 . 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 non-zero interaction effects. However, when running i-LDSC, only a percentage of the interacting SNPs {1%, 5%, 10%, 25% or 50%} were included in the estimation of . We once again generate traits with heritability such that the proportion of genetic variance explained by additive effects was equal to . 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 i-LDSC regression analyses. Overall, as discussed in the main text with results taken over 100 replicates, i-LDSC underestimated values of when there were unobserved interacting variants (see Figure 3—figure supplements 10 and 11). As expected, estimates of the additive variance component , 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 () and interaction () coefficients in the generative model of complex traits could potentially bias the additive and non-additive variance component estimates in LDSC and i-LDSC. 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 MAF-dependent procedure described above (i.e. ). We next selected a subset of the causal variants to be in cis-interactions. 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 (see Figure 3—figure supplement 12). In the second, we simply amplified the interaction effects to be a linear function (Figure 3—figure supplement 13A and C) or a squared function (Figure 3—figure supplement 13B and D) of the additive effects where . While testing 100 replicates for each value of , 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 broad-sense heritability levels = 0.3 and 0.6.
Preprocessing for the UK Biobank and BioBank Japan
Request a detailed protocolIn order to apply the i-LDSC 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 self-identified ancestral cohorts using data field 21000, unrelated individuals were selected by randomly choosing an individual from each pair of related individuals. This resulted in 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 Hardy-Weinberg equilibrium, using the PLINK --hwe 0.000001 flag to remove all SNPs with a Fisher’s exact test p-value . Finally, all SNPs with missingness greater than 1% were removed using the PLINK --mind 0.01 flag.
We then performed a genome-wide association study (GWAS) for each trait in the UK Biobank on the remaining 8,981,412 SNPs. SNP-level 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 --indep-pairwise 100 10 0.1. Using the parameters --indep-pairwise removes all SNPs that have a pairwise correlation above 0.1 within a 100 SNP window, then slides forward in increments of ten SNPs genome-wide.
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 i-LDSC scores from a reference panel. We pruned the reference panel using the PLINK command --indep-pairwise 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 i-LDSC GitHub repository (https://github.com/lcrawlab/i-LDSC). Using summary statistics from BioBank Japan, with scores calculated from the EAS population in the 1000 Genomes, we obtained i-LDSC heritability estimates for each of the 25 traits.
Data availability
Source code and tutorials for implementing interaction-LD score regression via the i-LDSC package are written in Python and are publicly available online on GitHub (copy archived at Crawford and Smith, 2024). Files of LD scores, cis-interaction 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 s-LDSC 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 (Bulik-Sullivan et al., 2020). When applying s-LDSC, 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 non-additive 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
-
SoftwareInteraction-LD score (I-LDSC) 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 genome-wide association studies and other large-scale problemsThe Annals of Applied Statistics 5:1780–1815.https://doi.org/10.1214/11-AOAS455
-
Estimation of non-additive 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 genome-wide 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/s41588-021-00827-w
-
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/s41467-020-17576-9
-
Quantifying the contribution of dominance deviation effects to complex trait variation in biobank-scale dataAmerican Journal of Human Genetics 108:799–808.https://doi.org/10.1016/j.ajhg.2021.03.018
-
PLINK: a tool set for whole-genome association and population-based 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/1297-9686-43-25
-
The influence of higher-order epistasis on biological fitness landscape topographyJournal of Statistical Physics 172:208–225.https://doi.org/10.1007/s10955-018-1975-3
-
Rare-variant 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 genome-wide 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 large-scale multiple regression with summary statistics from genome-wide association studiesThe Annals of Applied Statistics 11:1561–1592.https://doi.org/10.1214/17-aoas1046
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 (DMS-1439786)
- 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 (DBI-1452622)
- 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. DMS-1439786 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 DBI-1452622 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 EDC-induced epimutations are introduced or subsequently propagated, whether there are cell type-specific 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 cell-like (PGCLC) cells to BPS and found that differential incidences of BPS-induced 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 EDC-induced 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 sequence-to-function mapping remains challenging because neither function nor substrate specificity of biosynthesis enzymes can accurately be predicted. Here, we developed a knowledge-guided bioinformatic pipeline to solve these issues. We analyzed 1928 genomes of Pseudomonas bacteria and focused on iron-scavenging 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 iron-loaded 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.