Joint genetic analysis using variant sets reveals polygenic gene-context interactions.

Joint genetic models for multiple traits have helped to enhance association analyses. Most existing multi-trait models have been designed to increase power for detecting associations, whereas the analysis of interactions has received considerably less attention. Here, we propose iSet, a method based on linear mixed models to test for interactions between sets of variants and environmental states or other contexts. Our model generalizes previous interaction tests and in particular provides a test for local differences in the genetic architecture between contexts. We first use simulations to validate iSet before applying the model to the analysis of genotype-environment interactions in an eQTL study. Our model retrieves a larger number of interactions than alternative methods and reveals that up to 20% of cases show context-specific configurations of causal variants. Finally, we apply iSet to test for sub-group specific genetic effects in human lipid levels in a large human cohort, where we identify a gene-sex interaction for C-reactive protein that is missed by alternative methods.

Introduction

Understanding genetic interactions with external context (GxC), including environment, is a major challenge in quantitative genetics. Linear mixed models (LMMs) have emerged as the framework of choice for many genetic analyses, mainly because the random effect component in this class of models provides robust control for population structure [1, 2] and other confounding factors [3-5]. More recently, random-effect models have also been shown to be effective to test for polygenic effects from multiple causal variants that are in linkage [6-9] (variant sets). Additionally, multivariate formulations of LMMs have been developed to test for genetic effects across multiple correlated traits, both in single-variant analyses [10, 11] and more recently for joint tests using variant sets [12]. However, these existing multivariate LMMs have primarily been designed to increase the statistical power for detecting association signals, whereas methods to test for interactions are only beginning to emerge [10, 13].

Classical single-variant models for GxC use fixed effects to test for differential effect sizes of individual variants between contexts, either using an ANOVA [14-16] or LMMs [10, 17]. The main advantages of set-based tests compared to single-variant models are twofold. First, set tests reduce the effective number of tests and can account for effects due to multiple causal variants, thus increasing power for detecting polygenic effects [7, 8, 12, 18]. Second, we here show that joint tests across multiple contexts and sets of variants allow for characterizing the local architecture of polygenic-GxC interactions.

One way to generalize single-variant interaction tests to variant sets is using a model that assumes that context differences cause the same fold-differences in effect size across all genetic variants, such that all genetic effects in one context are proportional to the effects in a second context; a criterion that has also been considered to assess co-localization of multiple traits [19] (, middle). We denote this class of interactions rescaling-GxC. More generally, however, there may also be differences in the configuration of causal variants between contexts (, right), such that not all genetic variants show the same fold-difference between contexts, as some become more prominent in particular contexts and others less so. We denote these complex interactions heterogeneity-GxC. These two classes of interactions have different functional implications–the former suggest no difference in causal variants between contexts, and the latter suggest otherwise. Distinguishing between them is only possible using multi-variant models such as set tests, and is important for identifying different potential causal variants in different contexts.

Illustration of the iSet model and different architectures of genotype-context interactions.

(a) Alternative genetic architectures that are explicitly modeled in iSet: persistent effects, where causal variants have identical effects across contexts (left panel), rescaling-GxC effects, where the effects of causal variants in one context are proportional to those in a second contexts (middle), and heterogeneity-GxC effects, with changes of causal variants or their relative effect sizes between contexts (right). (b) Illustration of the multivariate linear mixed model (LMM) that underlies iSet. Model comparisons of LMMs with different trait-context covariance of the set component are used to define tests for general associations (mtSet), interactions (iSet) and heterogeneity-GxC effects (iSet-het). Additionally, the model can be used to estimate the proportion of variance that can be attributed to the corresponding genetic architectures (Methods). (c,d) Applications of iSet to a small simulated region. The total genetic effect was simulated as the sum of contributions from three loci with a persistent (left), rescaling-GxC (middle) and heterogeneity-GxC effects (right). (c) Manhattan plots of P values from a single-variant LMM [10] to test for associations (mtLMM) or interactions (mtLMM-int). Lower panel: Corresponding Manhattan plots for P values from set tests, considering a test for associations (mtSet), interactions (iSet) or heterogeneity-GxC (iSet-het), using consecutive regions (30 kb regions; step size 15 kb). Horizontal lines correspond to the α = 0.10 significance threshold (Bonferroni adjusted). P values of set tests are bounded (>10−6) by the number null model simulations to estimate significance levels (Methods). (d) Proportion of variance attributable to persistent effects, rescaling-GxC and heterogeneity-GxC, considering the same regions as in c.

We here propose a multivariate LMM to test for nteractions test between of genetic variants and categorical contexts (iSet) and to distinguish between rescaling-GxC and heterogeneity-GxC. We find that iSet yields increased power for identifying interactions and uniquely is able to robustly differentiate between rescaling-GxC and heterogeneity-GxC. We first validate iSet using simulations before applying the model to test for gene-by-sex interactions in blood lipid levels [20] as well as gene-by-environment interactions in an expression quantitative trait loci (eQTL) study [21]. We identify up to 20% of the stimulus-specific eQTLs as cases of heterogeneity-GxC, suggesting that context-specific causal variants are common.

Results

A mixed model approach to test for polygenic GxC

iSet generalizes previous multi-trait set tests [12], while considering the same trait measured in two (environmental) contexts. For a fully observed design, where the trait is measured in N individuals and each context, the phenotype matrix is modeled as the sum of a genetic effect from a set component and residual noise:

Here, and denote the design and the effect size matrices of additional fixed effect covariates and and are random effects that follow matrix-variate normal distributions: where corresponds to a local realized relatedness matrix [22] of the set of interest s, and denotes a diagonal covariance, which corresponds to independent and Identically distributed residuals. The trait-context covariance matrices and model correlations between contexts due to the set component (), and residual noise ().

A key insight derived here is that different assumptions on the structure of the trait-context covariance correspond to alternative genetic architectures that can be explained by a polygenic model (, Methods). Persistent genetic effects across contexts (no GxC) can be modeled using an LMM with a constant block covariance (); rescaling GxC, where effect sizes in different contexts are proportional to each other, can be captured by a trait-covariance with rank one. Note that genetic effects that act only in one context are a special case of this model and corresponds to a zero-rescaling coefficient. Finally, the most general architectures with different relative effect sizes between contexts (heterogeneity-GxC) can be captured by an LMM with a full-rank trait-context covariance (Methods). By comparing LMMs with these alternative covariance structures, it is possible to define set tests for general associations (mtSet), which identifies both persistent and context-specific effects, a test for genetic interactions, both with or without changes in the configuration of causal variants (iSet), and finally a test for heterogeneity-GxC effects (iSet-het), which is specific to differences between contexts that cannot be explained by rescaling ().

These multivariate LMMs can be fit using principles that were previously derived for multivariate set tests [12], and hence, provided the computations are suitable parallelized (Methods), iSet can be applied to large cohorts with tens of thousands individuals (). Permutation schemes are not well defined for interaction models [23], so we use a parametric bootstrap procedure [23] to estimate P values. An important advantage compared to previous interaction tests [13, 24–28] (Methods), is that iSet can be applied both to study designs where all individuals have been phenotyped in each context and when stratifying populations into distinct subgroups using a context variable (). iSet also provides control for population structure, either using principal components that are included as fixed covariates, or using an additional random effect term (Methods). Finally, iSet can also be used to estimate the total phenotypic variance explained by variant sets and the relative proportions captured by persistent, rescaling-GxC and heterogeneity-GxC effects (Methods).

To illustrate the polygenic interactions that can be detected using iSet, we first considered a basic simulated example (). We simulated genetic effects for one quantitative trait in two contexts, considering polygenic effects at three distinct loci (Methods): a region with persistent genetic effects, a region with rescaling-GxC and a region with heterogeneity-GxC effects. We tested consecutive regions (30kb region, 15kb step) using the three tests provided by our model (mtSet, iSet, iSet-het), finding that by combining these results, it was indeed possible to resolve the architecture of each of the simulated regions (). In particular, this example illustrates that, unlike single-variant tests, iSet-het can be used to discern heterogeneity-GxC effects specifically.

Simulated data

Next, we used simulations based on genotypes from the 1000 Genomes project [29] to assess the statistical calibration and power of iSet. We generated a population of 1,000 individuals based on genotype data from European populations, initially simulating one quantitative trait measured in two distinct contexts in all individuals (Methods).

First, we considered data with simulated persistent polygenic effects, confirming that both iSet and iSet-het are calibrated when no interaction effects are simulated (, ). Analogously, we also confirmed that iSet-het is calibrated when only rescaling-GxC effects are considered (), and we assessed the robustness of iSet to different types of model misspecification ().

Simulated data to assess statistical calibration and power of iSet.

(a) QQ plot for the P values obtained from iSet and iSet-het when only persistent genetic effects were simulated. The step in the QQ-plot for large p-values is observed because the trait-context covariances are required to be positive-semidefinite. (b) Power comparison of alternative models for detecting simulated interactions, considering rescaling-GxC effects (without heterogeneity-GxC) for increasing numbers of simulated causal variants at constant total genetic variance. Compared were iSet and a single-variant interaction test (mtLMM-int) [10], using two alternative approaches to adjust for multiple testing of single variant methods (Bonferroni or eigenMT). (c) Lower panel: analogous power comparison as in b, when varying the proportionality factor of effect sizes between contexts. A proportionality factor of zero corresponds to genetic effects that act only in one of the contexts. See for the relationship of the proportionality factor and fold differences. iSet-het was not considered, because all simulated rescaling-GxC are consistent with the null model of iSet-het. Top panel: average fraction of genetic variance attributable to persistent, rescaling-GxC and heterogeneity-GxC effects for the corresponding simulations. (d) Analogous comparison as in c but for simulated heterogeneity-GxC effects, when varying the correlation of the total genetic effect between contexts. Additionally, we also considered iSet-het to test for heterogeneity-GxC, which was best powered for heterogeneity-GxC effects that were uncorrelated between contexts. White stars denote default parameter values that were kept constant when varying other parameters (). Statistical power was assessed at 5% FDR across 1,000 repeat experiments.

We compared iSet to single-variant interaction tests [10] (mtLMM-int) (Methods), considering a wide range of different settings ( Methods). Because single-variant methods perform one test for each variant in the set (), we adjusted for multiple testing using one of two approaches: i) conservative Bonferroni adjustment (Bonferroni) or ii) a recently proposed method that estimates the effective number of independent tests based on the local structure of linkage disequilibrium (LD) [30] (eigenMT). Note that existing set-based interaction tests cannot be applied to complete designs with repeat measurements and hence were not considered (Methods, ); see below and for additional experiments where these methods were used. As expected, the power advantages of iSet compared to single-variant models were largest when multiple causal variants were simulated (, for constant total genetic variance, Methods). However, iSet was better powered than mtLMM-int even for a single causal variant. Identical simulations based on synthetic independent genotypes () revealed that this effect is predominantly due to local LD and advantages due the reduced number of total tests. We also considered the impact of different proportionality factors of genetic effects between contexts. All models were best powered to detect GxC for negative proportionality factors (opposite effects), or when the proportionality factor was close to zero (context-specific effects) ().

Next, we simulated traits with context-specific causal variants (heterogeneity-GxC). Heterogeneity-GxC is detectable when there is a change in the local causal configuration, which corresponds to the absolute correlation of local genetic effects between contexts (r) smaller than 1; the greater the heterogeneity GxC effects, the smaller the absolute correlation. Presence of GxC effects under tightly correlated genetic effects (r ≈ ±1) cannot be distinguished from rescaling-GxC. To simulate these different settings, we randomly selected two causal variants in each context and varied the extent of correlations of the genetic effect between contexts (). When using iSet-het for detecting heterogeneity-GxC effects, the model was best powered when there is a moderate to large change in causal configuration, corresponding to low correlated genetic effects (>70% power for r2 < 0.16, ). We also considered additional settings with larger numbers of causal variants (), and we assessed the accuracy of iSet-het to classify interaction effects into heterogeneity-GxC or rescaling-GxC effects (, Methods). Taken together, these results confirm that iSet-het is a robust test for heterogeneity-GxC.

We also investigated the proportion of local genetic variance that can be explained by models with persistent, rescaling-GxC and heterogeneity-GxC for the corresponding simulations (). The persistent effect model explained large proportions of the simulated genetic variance, even in the presence of positively correlated GxC, but could not capture variance due to GxC effects with negative rescaling (). An LMM that models rescaling-GxC did account for negative and positive rescaling, and captured some of the heterogeneity-GxC effects (). Finally, variance contributions that were exclusively captured by a heterogeneity-GxC model were largest for uncorrelated context-specific genetic effects, the same regime where the corresponding test is best powered (). We also confirmed that the most flexible heterogeneity-GxC model yields unbiased estimates of the total genetic variance in genomic regions, whereas other models were biased for some simulated architectures ().

Finally, we considered simulations where we varied both the size of the testing region and the simulated causal region, using a sliding window analysis (, Methods). Overall, iSet was markedly robust to the window size, and was best powered when the sizes of the testing region approached the size of the simulated causal region, which is in line with previous findings for set-based association testing [12]. We also observed that iSet-het is best powered for small causal regions (up to 100kb), and the power for detecting heterogeneity-GxC deteriorated when analyzing larger regions.

Analysis of stimulus-specific eQTLs in monocytes

We next applied iSet to test for stimulus-specific genetic effects in a monocyte stimulus eQTL study [21]. We considered gene expression profiles for 228 individuals in four stimulus contexts: naive state (no stimulation), stimulation with interferon-γ for 24 hours (IFN), and stimulation with lipopolysaccharide (LPS) for two and 24 hours.

We applied iSet to test for pairwise interaction effects, considering the naive monocyte state and each stimulus condition in turn, performing a single test using proximal cis acting variants (plus or minus 50kb from the transcription start site; Methods). After quality control, we considered 12,677 probes and tested for cis associations (mtSet), GxC interactions (iSet) and for heterogeneity-GxC effects (iSet-het). For comparison, we also considered a conventional multi-trait LMM [10] and tested for associations and interactions in the same genomic regions, using eigenMT [30] to adjust for multiple testing (Methods). Although there was substantial overlap of the probes and stimulus conditions for which different methods identified significant interactions (), iSet was better powered (32.7% power increase; 5,068 versus 3,818 probes and stimuli with an interaction; FDR<5%, ). Additionally, iSet-het identified 1,135 probes and stimulus contexts with significant heterogeneity-GxC effects (). This shows that a substantial proportion of stimulus-specific eQTLs are associated differences in the configuration of causal variants, suggesting context-specific regulatory architectures.

Analysis of stimulus-specific eQTLs in monocytes.

(a) Number of probes with at least one significant cis association (Association test) or genotype-stimulus interaction (Interaction test) for alternative methods and stimulus contexts. Considered were the proposed set tests (mtSet, iSet, iSet-het) as well as single-variant multi-trait LMMs (mtLMM, mtLMM-int [10]), testing for genetic effects in cis (100kb region centered on the transcription start site; FDR < 5%). Additionally, iSet-het was used to test for heterogeneity-GxC effects. Individual rows correspond to different stimulus contexts with “All” denoting the total number of significant effects across all stimulus contexts. (b) Venn diagram of probes and stimuli with significant interactions identified by alternative methods and tests (across all stimuli). (c) Bivariate plot of the variance attributed to persistent genetic effects versus genotype-stimulus interactions for all probes and stimuli. Significant interactions are shown in red. Density plots along the axes indicate the marginal distributions of persistent genetic variance (top) and variance due to interaction effects (right), either considering all (black) or probe/stimulus pairs with significant interactions (iSet in a, dark red). (d) Average proportions of cis genetic variance attributable to persistent effects, rescaling effects and heterogeneity-GxC, considering probe/stimulus pairs with significant cis effects (5% FDR, mtSet), stratified by increasing fractions of the total cis genetic variance. Shown on top of each bar is the number of instances in each variance bin. The top panel shows the density of probes as a function of the total cis genetic variance.

Although on average the proportion of variance explained by GxC tended to be smaller than for persistent effects (median 3.7% for GxC versus median 9.5% for persistent effects, for probes with significant GxC, ), GxC was the driving genetic source of variation for 11.8% of the significant cis eQTLs (; defined as explaining 50% or more of the cis genetic variance). Consistent with previous reports [31, 32], we observed that genes with large relative GxC effects were associated with weak overall cis effects, whereas eQTLs with large effect sizes tended to be persistent across stimuli ().

Mechanistic underpinning of heterogeneity eQTLs

To better understand the mechanisms that underlie genes with detected heterogeneity-GxC effects, we used an LMM with step-wise selection [33], identifying 15,756, 2,690 and 457 eQTLs (across all probes and contexts) with a single significant association, significant secondary and significant tertiary associations respectively (FDR < 5%, Methods, ). Probes with significant heterogeneity-GxC were more likely to harbor multiple independent associations (), confirming that heterogeneity-GxC eQTLs have complex genetic architectures.

Characterization of genes with significant heterogeneity GxC for stimulus eQTLs in monocytes.

(a) Cumulative fraction of probe/stimulus pairs with increasing numbers of distinct univariate eQTLs (average of the naïve and the stimulated state using step-wise selection) for different gene sets (Methods). Shown are cumulative fractions of all probe/stimulus pairs (All), those with significant cis associations (mtSet), pairs with significant GxC (iSet) and instances with significant heterogeneity GxC (iSet-het). (b) Breakdown of 1,281 probe/stimulus pairs with significant heterogeneity GxC into distinct classes defined using the results of a single-variant step-wise LMM (Methods). (c-e) Manhattan plots for representative probes with significant heterogeneity GxC effects. Grey boxes indicate the gene body. (c) Manhattan plot (left) and χ2 statistics for variants in both contexts (right) for the gene SLC1A4. Dark circles indicate distinct lead variants in both contexts (r2<0.2). (d) Manhattan plot after conditioning on the lead variant (secondary associations in the stepwise LMM) for the gene PROK2. The star symbol indicates the shared lead variant in both contexts. The conditional analysis reveals a secondary association that is specific to the naïve state. (e) Analogous plot as in c for the gene NSUN2, for which the single-variant model did not provide an interpretation of heterogeneity-GxC. (f) Breakdown of probe / stimulus pairs with shared lead variants, stratified by concordance of the effect direction (opposite-direction versus same-direction eQTLs) and significance of the heterogeneity-GxC test (heter vs No heter; FDR 5%). eQTLs with opposite effects were enriched for significant heterogeneity-GxC (2.2 fold enrichment, P<4e-2).

When overlaying heterogeneity GxC eQTLs detected using iSet-het with the results obtained from the single-variant step-wise LMM, we could attribute 46.2% of the heterogeneity-GxC effects (524 out of 1,135) to context-specific lead variants (defined using r2<0.2, FDR<5%, , Methods). For an additional 14.6% of the heterogeneity eQTLs (166/1,135) the lead variants from a single-variant analysis were in high LD (r2>0.8), with context-specific secondary effects (). The remaining 445 heterogeneity eQTLs (39.2%) could not be annotated using single-variant models.

One reason why heterogeneity GxC effects cannot be annotated using a single-variant model are differences in power. Indeed, for 22.2% of the heterogeneity-GxC cases without a single-variant interpretation (99/445), the single-variant LMM did not yield a significant effect in either of the two contexts (). For an additional 58.2% of the unannotated heterogeneity GxC effects (259/445), the single-variant LMM lead variants were in weak linkage (0.2

Finally, we explored the relationship between probes with heterogeneity GxC and opposite effects as defined using conventional single-variant models. We classified associations as opposite effects when context-specific lead variant were in high LD (r2>0.8) and the effect on gene expression was in different directions (Methods). This approach identified 67 eQTLs with reversed effect directions between contexts. iSet-het detected significant heterogeneity-GxC for 8 of these eQTLs, a 2.2 fold enrichment (P<5e-2) compared to eQTLs with consistent effect directions between contexts (238 gene/stimulus pairs with significant heterogeneity-GxC out of 4,119 eQTLs with consistent direction, ). Similar enrichments were also observed when considering individual stimulus contexts, resulting in significant enrichments for two out of three stimulus contexts (P<5e-2, fold change>4 in naïve/IFN and naïve/LPS-24h, ). Among the genes with significant heterogeneity-GxC are OAS1, LMNA and PTK2B, opposite-effect eQTLs that have been reported in the primary analysis of the same data [21] ().

Using iSet to test for interaction effects in stratified populations

Thus far, we have considered settings with repeat measurements, where the same phenotype is measured in all individuals and contexts. Next, we considered applications of iSet to studies where individuals are phenotyped in only one of the two contexts (, Methods). This is a common strategy in investigation of genotype-context interactions, where a population is stratified using a context variable.

We considered simulations analogous to those for complete designs () to validate iSet for this design. We again confirmed statistical calibration of iSet () and found similar power benefits as for complete designs (, ). In addition to single-variant LMMs, we also compared to a recently proposed set test for interactions (GESAT; [13]), which is designed for stratified populations. Notably, iSet was consistently better powered than GESAT, most likely because GESAT does not model correlations of the local genetic effect between contexts (Methods, ).

Application of iSet to stratified designs.

(a,b) Power comparison of iSet and alternative methods using simulated data where each individual is phenotyped in one of two contexts. Shown is a comparison of power for alternative methods. (a) Power to detect interactions when simulating rescaling-GxC for increasing numbers of causal variants. (b) Power when varying the factor of proportionality of the variant effect sizes between contexts. Considered were iSet, a single-variant interaction test (mtLMM-int, [10]) as well as the interaction sequence kernel association test (GESAT, [13]), a set test designed for stratified populations. For single-variant models, two alternative approaches to adjust for multiple testing were considered (Bonferroni, eigenMT). (c) QQ-plot of P values from genotype-sex interaction tests for C-reactive protein levels using individuals from the Northern Finland Birth Cohort [20], considering the same methods.

Next, we applied iSet to test for genotype-sex interactions in four lipid-related traits (fasting HDL and LDL cholesterol levels, triglycerides and C-reactive protein) measured in 5,256 unrelated individuals from the Northern Finland Birth Cohort (NFBC1966 [20]). We tested consecutive 100kb regions (step size 50 kb; 52,819 genome-wide tests), and compared iSet to GESAT and the single-variant interaction test (Methods).

iSet retrieved one genome-wide significant interaction (C Reactive protein, chr1:40,450,000; P = 1.47x10-6; FWER<10%), whereas alternative set tests and the single-variant models did not yield significant effects (, ), even when using dense genotypes derived using imputation strategies (). This interaction was located within 400kb of an interaction identified in a large meta study (66,185 individuals [34]), which reports both an association for C-reactive protein at the same locus (P<6x10-11) as well as a nominally significant interaction with sex (P<5x10-3). Finally, a local single-variant analysis, separately for female and male individuals, provided evidence that this interaction reflects a male-specific genetic effect ().

iSet revealed a second suggestive interaction with sex for LDL cholesterol levels (chr3:121,850,000, ). Although this effect failed genome-wide significance (FWER<20%), iSet again yielded stronger evidence than other methods (PiSet = 3.7x10-6, PGESAT = 4.8x10-6, PmtLMM-int = 3.2x10-5). Among the genes at this locus is ADCY5, which has been linked to blood glucose levels in large meta-analyses [35, 36] and hence is a plausible candidate to affect LDL via glucose regulation [37].

Finally, we note that context stratification of quantitative traits can increase power for detecting associations rather than interactions, which is similar to previous strategies applied for single-variant analyses of quantitative [38] and categorical traits [39, 40]. Using this generalized association test, we identified three additional associations that were missed by conventional set tests and other methods (). These include the same locus with a sex-specific effect on C-reactive protein (chr1:40,450,000, P = 1.42x10-7 using mtSet, P = 1.89 x10-3 using a standard set test), and two associations for HDL cholesterol levels and triglycerides, both of which were replicated in larger meta analyses [41].

Discussion

We have here proposed iSet, a method based on linear mixed models to test for gene-context interactions using variant sets. On simulated data as well as in applications to gene expression and human lipid-related traits, we have demonstrated that iSet yields increased power and improved interpretation for interaction effects compared to previous methods.

Methods for the joint analysis of multiple traits, including tests for genetic interactions, are not new per se. Most previous studies have used set-based methods to test for associations [7, 8, 12, 18], whereas tests for genotype-context interactions are still primarily carried out using single-variant models [10, 17]. iSet unifies several previous models (Methods), and uniquely offers set-based interaction tests on phenotypes in different contexts under the same or different (stratified) genetic backgrounds. Additionally, we have shown that set-based interaction tests can be useful to disentangle the genetic architecture of such loci, discerning consistent changes of genetic effects between contexts (rescaling-GxC) and changes in the configuration of causal variants (heterogeneity-GxC). The heterogeneity GxC test we propose is related to co-localization tests [19, 42, 43], however with a different objective.

In applications to a stimulus eQTL study, iSet has yielded increased power compared to single-variant tests (), and we have shown that approximately 20% of the gene-stimulus interactions are associated with significant heterogeneity-GxC. This suggests that changes in the genetic architecture between stimulus contexts are relatively common. Additionally, we have observed that genes with opposite effects are enriched for heterogeneity-GxC. This finding points to a possible bias whereby opposite effects identified using single-variant models may in part be due to context-specific causal variants that are LD-tagged by a shared lead variant. Notably, although iSet is better powered overall, there may be true interactions that can be detected using single-variant models and are missed by set tests (). Hence, iSet should be considered as a complementary method and not as a replacement of conventional single-variant tests.

The proposed iSet model is not free of limitations. First, scalable inference in our model is achieved by exploiting the low-rank structure of variant sets, meaning that the number of variants in the analyzed region is typically small compared to the number of individuals. Similar to previous set-based tests [12], there are trade-offs between power and resolution, in particular when analyzing data from densely imputed or sequenced cohorts. General strategies for the design of optimal testing regions, for example using genome annotations and LD information, are an important area of future work. iSet is computationally efficient in cohorts with fully observed designs, or when stratifying a cohort using a context variable. Intermediate designs, for example in fully observed designs with missing phenotypes, may also be considered, however currently require the use of separate imputation schemes [11, 44]. It is also worth noting that the test for heterogeneity-GxC (iSet-het) will be most accurate if all individuals are phenotyped in each context. Although in principle the model can also be used in stratified designs, there may be concerns that false positive heterogeneity GxC effects can arise due to technical factors, for example due to differences in genotyping accuracy or variant allele frequencies in the corresponding sub populations. A related issue is the need to choose the size of the region-set appropriately. While we find that the model is overall robust across a wide range of region sizes (), the model will be best powered if the size of true causal regions approximately matches the testing region size, in particular for identifying heterogeneity-GxC effects.

Finally, we have here focused on pairwise analyses of different contexts. In principle, the model could also be applied to analyze multiple related context and different traits, and the model could be extended to handle continuous environmental states, which currently require discretization. A related extension of the model is to test for genetic effects that are exclusive to one of the considered contexts. Developments in these directions are future work.

Methods

Software availability

iSet is freely available as part of the LIMIX package (https://github.com/limix/limix). Tutorials for using iSet either as command line tool or via a Python API can be found at https://github.com/limix/limix-tutorials/tree/master/iSet.

The interaction set test (iSet)

To derive the model, we start assuming a fully observed design, where phenotypic measurements are available for all individuals and for each context. Briefly, the N × C phenotype matrix Y for N individuals and two or more contexts (C) is modeled as sum of fixed effects of K covariates, effects from S genetic variants in the region of interest (set component) and residual noise:

Here (N × K) and (N × S) denote respectively the fixed-effect covariates and the standardized genotypes of the variant set and (K × C) and (S × C) denote the corresponding effect sizes. The noise component is assumed to follow a matrix-variate normal distribution, ∼MVN(0,,), where is a C × C covariance matrix that models residual covariances between traits. Note that in this formulation, population structure can be accounted for by including the leading principal component of the N × N (global) realized relatedness matrix [22] into the model as fixed effects [12]. In human populations, 10–20 principal components are typically sufficient to adjust for such structure [45]. Note that iSet can also account for population structure using an additional random effect term into the model (see ). While computationally more expensive, this approach provides for additional robustness and calibration when analyzing cohorts with related individuals (see [12] for a discussion). All experiments reported here have been carried out using adjustment based on principal components, considering 10 PCs.

Relationship between the trait-context covariance and the genetic architecture

To simplify the notation, we consider the case of two contexts, however the same derivation holds for larger numbers of contexts. Different genetic architectures between contexts are cast as specific assumptions on the S × C matrix of the variant effect sizes . A persistent genetic effect can be expressed as , where a shared genetic signal (S × 1) has a common scale a0 in both contexts. Rescaling-GxC can be expressed as = , where a common genetic signal (S × 1) is modulated by context-specific scales = [a1; a2] (2 × 1). Finally, in the most general case, the configuration of causal variants is independent between contexts, corresponding to = [1, 2], with genetic signals 1 (S × 1) and 2 (S × 1) and general scaling factor matrix = [a11, a12; a21, a22] (2 × 2).

Marginalizing over the genetic signal , 1 and 2, assuming independent unit-variance normal prior distributions, results in a marginal likelihood of the form

Here, vec denotes the stack-column operation, ⊗ the Kronecker product, is the C × C trait-context covariance for the set component and is the local realized relatedness matrix (). The stated alternative generative models for the structure of have a one-to-one correspondence with alternative covariance structures for (), where for persistent effects is a block covariance ( = a012×2), rescaling-GxC correspond to with a rank-one structure ( = ) and the most case of independent architectures in both contexts can be captured by a full rank covariance ( = ). In order to test for associations, we additionally consider a null model without a set component ( = 0). The trait-context covariances , and the fixed effect weights are estimated using (restricted) maximum likelihood with constraints to obey the alternative structures of . Model parameters are optimized based on the restricted log marginal likelihood as objective, using a low-memory Broyden-Fletcher-Goldfarb-Shanno optizmier (L-BFGS) [46], implemented in the fmin_l_bfgs_b optimisation method of the SciPy python library. Specific tests are implemented as pairwise likelihood ratio (LLR) tests, considering LMMs with different trait-context covariances (See also ):

mtSet: full-rank versus null

iSet: full-rank versus block covariance

iSet-het: full-rank versus rank-one

iSet support these tests in fully observed designs and in stratified populations. Care is required when applying the heterogeneity-GxC test to stratified designs. Differences in allele frequencies between strata and other confounding factors could potentially result in false positive heterogeneity-GxC signals.

Obtaining P values

Empirical P values are estimated from the distribution of LLRs under the null. As permutation procedures are not well defined for interaction tests, for both iSet and iSet-het we generate test statistics from an empirical null distribution using a parametric bootstrap procedure [23]. Briefly, this procedure consists of sampling phenotypes from the null model with parameter values that maximize the likelihood on the observed data. Similarly to [8, 12], we consider a small number of parametric bootstraps for each region (typically 10–100 bootstraps) and pool the obtained null LLRs across all tested regions. The estimated distribution of null LLRs is used to obtain empirical P values. In an analysis of T genomic regions, the procedure to obtain P-values for the three tests can be summarized as follows:

fit the no-association model (null), the block covariance model (block), the rank-one covariance model (rank-one) and the full-rank covariance model (full) and estimate LLRs for mtSet (full vs null), iSet (full vs block) and iSet-het (full vs rank-one);

for each genomic region, sample J LLRs from the null for each of the three tests (J permutations for mtSet, J parametric bootstraps for iSet and J parametric bootstraps for iSet-het);

for each of the these tests, pool the JT null LLRs across regions to obtain an empirical null and compute empirical P values.

Note that the number of parametric bootstraps/permutations will determine the minimum P value that can be obtained. For example, for T tests and J = 30 bootstraps the minimum P value that can be estimated is 1/(JT), which correspond to a FWER of 1/J ≈ 0.03. While 30–100 bootstraps will be sufficient to reach typical thresholds in genome-wide studies, more stringent thresholds on significance levels (FWER< = 1%) require a larger number of parametric bootstraps (see section below for computational considerations). For mtSet we use the same procedure but with permutations [12].

Data design, relatedness and scalability

Parameter inference using naïve implementations to fit the marginal likelihood model in iSet (Eq (2)) would scale cubically with the number of samples and contexts. iSet is optimized for cohorts with unrelated individuals, in which case population structure can be accounted for by including the top principal components (PCs) as fixed effect covariates. We have adapted prior work to multi-trait set tests, which for fully observed designs results in a computational complexity of , where N denotes the number of individuals, R denotes the number of variants in the region, NPC is the number of PCs and t corresponds to the number of function evaluations of the optimizer (See [12] for details). iSet will be most efficient when the number of variants in the set is small compared to the number of individuals. To enable applications to stratified cohorts, we have extended this inference scheme to designs where phenotype data from each sample are observed in only one of the contexts, resulting in a computational complexity of (see ). For cohorts with related individuals and fully observed designs, iSet can also be applied with an additional random effect term in the model. In this case, we again re-use efficient inference schemes for multi-trait set tests [12]. This model is computationally more expensive, and has computational complexity O(N3 + N2R + tNR2). Notably, the computational scaling as stated assumes that the number of variants in the testing set R is constant, and in particular does not increase as a function of the sample size N. In practice, however, there may be an implicit dependency between N and the number of variants R, for example when low frequency and rare variants are included in the analysis. Consequently, the empirical compute cost as a function of N could be larger than linear. A second setting are small cohorts, where the number of variants R could exceed N. In order to retain efficiency in this setting, iSet detects when the low rank assumption of the region set is violated and reverts to an eigen decomposition in the space of individuals, and hence the computational complexity is bounded by O(). See for a summary of the computational complexity of iSet for alternative data designs, strategies to adjust for confounding and analysis settings.

Moreover, although iSet is primarily designed for the analysis of pairs of contexts, the model can also be applied to multiple contexts. The computational cost of iSet as a function of the number of contexts C is O(tNC2 + tC3).

Empirical runtime estimates in were reported for different designs, using synthetic cohorts generated using data from the 1000 Genomes Project (phase 1, ). We report the average per-region compute time measured on 100 regions with size 30 kb, considering a single core of an Intel Xeon CPU E5-2670 2.60-GHz to fit iSet. The runtime for all three considered tests, mtSet, iSet and iSet-het (including bootstraps) on the eQTL analysis took on average 28.7s per gene, resulting ~100h of compute time for a genome-wide analysis using a single core. Similarly, the runtime for all three tests for the NFBC data was on average 107s per testing region, resulting in ~1,500h of compute time for a genome-wide analysis.

Finally, we note that while the linear scaling allows for application to larger datasets, the computational cost of iSet is typically larger than for a conventional single-variant model. Consequently, for genome-wide analyses it will be required to parallelize the computational operations across multiple compute cores. The software implementation of iSet supports the parallelization across multiple compute nodes and cores.

Variance decomposition model

The LMMs in iSet can also be used to estimate the phenotypic variance explained by the variant set for their persistent, rescaling-GxC and heterogeneity-GxC effects ().

Relationship to existing methods

iSet generalizes previous interaction set tests and multi-trait mixed models. Existing interaction set tests [13, 24–28] are designed for the analysis of stratified individuals and are not applicable to designs with repeat measurements, where the same trait is phenotyped in the same individuals in multiple contexts. Moreover, these existing methods do not account for correlated genetic effects within the region set and their underlying LMMs assume that the signal to noise ratio is identical in both contexts. The iSet model is more flexible and accounts for arbitrary genetic correlations and residuals covariances, using a null model that is similar to previous single-variant interaction tests [10]. iSet combines the advantages of several of these previous models; see for details.

Choice of the window size

As for any set test, the size of the region set is an important parameter in iSet. The specific choice will depend on the biological application, LD and marker density. We have previously explored trade-offs between the computational efficiency and power of association tests for different choices of the window size [12]. We here examined how the choice of the window size affects the power of detecting interaction and heterogeneity-GxC by considering sliding-window experiments with alternative windows sizes in simulations ( see below). For the simulation experiments shown in , we considered sets with 30kb, which captures local LD in the data (). For the analysis of the stimulus eQTL study, we have considered gene-based sets using a 100kb cis genetic region centered on the TSS, which is in line with other cis eQTL analyses [32]. Finally, for the genotype-sex interaction analysis in blood lipid levels we followed [12] and considered a sliding window approach with 100kb regions and a step size of 50kb.

Simulation study for fully observed designs

Simulations were carried out using a synthetic cohort of 1,000 individuals derived from genotypes of European populations in the 1000 Genomes project [29] (phase 1, 1,092 individuals, 379 Europeans). Following [12, 47], we composed synthetic genotypes as a mosaic of real genotypes from individuals of European ancestry, while preserving population structure (). We considered single-nucleotide polymorphism with a minor allele frequency of at least 2% (). In all simulations, we simulated two contexts, modeled as the sum of a genetic contribution from a 30kb causal region, effects due to population structure, hidden covariates and identically distributed Gaussian noise. Effects due to population structure and hidden confounders were simulated with partial correlations across contexts, explaining variable proportions of the total phenotypic variance in each context (, ).

Statistical calibration

To assess the calibration of P values obtained from the interaction test (iSet) and the test for heterogeneity-GxC (iSet-het), we considered 100,000 datasets with two contexts where only persistent genetic effects (no interactions) were simulated (). For each simulation we randomly selected a 30kb region and generated phenotypes simulating persistent effects from four causal variants and tested for GxC interaction in the region. To estimate P values, we used 30 parametric bootstraps for each test, resulting in a total of 3,000,000 null LLRs to estimate P values. Analogously, we assessed the calibration of iSet-het, where exclusively rescaling-GxC effects were simulated (). Again, we considered 30 parametric bootstraps for each test and pooled LLRs to estimate P values. We also assessed the calibration of iSet and iSet-het when simulating violations from the model assumptions (, ). These included (i) deviations from the polygenic model, using causal regions with a single causal variant, (ii) effects due to epistatic interactions of pairs of variants in the region and (iii) deviation from the Gaussian residual assumption by introducing outlying samples. Outliers were simulated by adding a heavy tailed noise contribution to 0.1% of the samples (effect ±8, explaining ~5% of the total sample variance on average). For simulated outliers, we compared iSet on un-normalized phenotype data as well when quantile normalizing the phenotypes to a normal distribution. Quantile normalization has previously been shown to be effective in combination with multi-trait LMMs [11, 12], and is the recommended default approach for iSet. Note that while quantile normalization does help to increase the robustness of multi-trait LMMs, it should be applied with caution when the noise is strongly heavy tailed or heteroskedastic. In such scenarios, quantile normalization could introduce false positive GxC signals by altering the relative scale of the traits.

Comparison with alternative methods

For comparison, we considered single-variant interaction tests as in [10] (mtLMM-SV-int), using an implementation in LIMIX [48]. To obtain region-based P values, we considered the minimum P value across all variants in the region, following adjustment for multiple testing. We consider two alternative strategies to adjust for multiple testing within variant sets: i) a conservative Bonferroni approach and ii) the recently proposed eigenMT model [30], which adjusts for the effective number of independent tests estimated based on the local LD structure. Existing set tests are not applicable for fully observed designs and hence were not considered ().

Power comparison

To assess power of iSet for alternative genetic architectures, we simulated interaction effects from a 30kb region either considering rescaling-GxC effects or more general effects that include heterogeneity-GxC, using the simulation settings in . The total variance explained by the causal region across all traits was set to 2%. In the case of rescaling-GxC, we varied i) the number of causal variants in the region (from 1 to 20; ), and ii) the proportionality factor of the effect sizes between contexts (from -1 to 1, ). When simulating general effects that include heterogeneity-GxC, we randomly selected an equal number of context-specific causal variants and monitored the correlation of the total simulated genetic effects across contexts, thereby controlling the extent of heterogeneity-GxC. Again, local genetic effects were simulated to explain 2% of the total phenotypic variance in each context. We varied (see ) i) the extent of simulated heterogeneity-GxC () and ii) the total number of causal variants across contexts (). For each parameter setting, we considered 1,000 repeat experiments. To obtain P values for set tests we considered 30 parametric bootstraps for each test and computed empirical P values from 30,000 null LLRs in each simulated scenarios. We used the Benjamini-Hochberg procedure to adjust for multiple testing across repeat experiments and assessed all methods in terms of power at a fixed FDR<5%.

Illustration case

For the simulated example region to illustrate iSet and alternative genetic architectures (), we used a simulation procedure analogous to the strategy described above. Phenotypes were simulated as the sum of genetic effects from three distinct causal regions (30kb) within a 5Mb region on chromosome 13, harboring respectively persistent, rescaling-GxC and heterogeneity-GxC effects. The effects from individual regions was simulated to explain 5% of the total phenotype variance.

Comparison of iSet-het with a baseline test for heterogeneity-GxC

As an additional assessment of the accuracy of iSet-het to detect heterogeneity-GxC effects, we tested how well the model discriminates between regions with and without simulated heterogeneity-GxC. We considered the identical 10,000 regions in for which no heterogeneity-GxC effects were simulated as well as the 10,000 regions in with heterogeneity-GxC. We ranked all 20,000 regions based on the LLR of the heterogeneity test and used the receiver-operating characteristic (ROC) and precision-recall curves () to assess the ability of discriminating between these types of genetic effects. For comparison, we also considered a univariate baseline approach, scoring regions with significant associations using the squared Pearson correlation between the lead variants in both contexts (low squared Pearson correspond to high rank). We considered alternative significance thresholds on region-based P values obtained using eigenMT (P < 0.5, 0.01, 1e-2, 1e-3).

Sliding window simulation experiments

To study the effect of alternative sizes of the testing region on the power of iSet and iSet-het under different simulated scenarios, we considered sliding-window experiments using simulated data, when varying both the size of the simulated causal and of the testing region. Phenotypes were simulated across two contexts using the approach as described above, considering a causal region with variable size (30kb, 100kb, 300kb and 1Mb). For each of these simulation settings, we carried out a sliding window analysis in the surrounding 1 Mb region with testing windows of 30kb, 100kb, 300kb and 1Mb (the step size was set to the half of the size of the testing region). We used Bonferroni to adjust for multiple testing across regions. For comparison, we also considered the single-variant test for interactions (mtLMM-int), applied to the same variants in the 1Mb region, and used eigenMT to adjust for multiple testing across variants while accounting for LD. For each scenario, we considered 200 repeat experiments and assessed power at FDR = 5%. We considered either simulated pure rescaling-GxC () or more general effects (rescaling+heterogeneity-GxC, ). For both sets of simulations, we considered the default simulation parameter values ().

Monocyte eQTL dataset

Data pre-processing

The dataset consists of gene expression levels from primary monocytes, both in a naïve state and three different stimulus contexts, profiled in 432 genotyped individuals of European ancestry. Gene expression levels in the naïve state, after exposure to IFN-γ, after 24-hour LPS and after 2-hour LPS were available for 414, 367, 322 and 261 individuals respectively. Normalization, correction for batch and probe filtering were done as in [21]. Following [21], we only considered probes that (i) map to only one genomic location, (ii) do not overlap with SNPs (MAF>1% in Europeans populations of 1000 Genomes Project), (iii) map to regions on autosomal chromosomes, and (iv) were detected in sufficient number of samples (see [21] for more details). Additionally, we discarded probes that could not be mapped to Ensembl gene IDs. Collectively, these filters resulted in 12,677 probes for analysis (out of 15,421). We further limited our analysis to the set of 228 individuals for which gene expression levels were available in all the four (stimulus) contexts. To account for hidden covariates and confounding factors, we applied PEER [49] with default parameter values, fitting 30 hidden factors across all samples (individuals and stimulus states). PEER residuals for each gene and context were quantile-normalized to a standard normal distribution and used for all genetic analysis. Again, following the primary analysis [21], genotypes were imputed against the 1000 Genomes Project reference panel. After excluding variants with MAF<4%, variants with low imputation score (<0.9) and variants that deviate from the Hardy-Weinberg equilibrium (pv<10−3), we were left with 5,729,118 genome-wide variants (4,967,901 unique variants).

eQTL mapping

Association and interaction tests were carried out considering 100 kb regions centered on the transcription start site of genes corresponding to individual probes (). All tests were applied considering a pair-wise approach, jointly testing for eQTLs in the naive state and one of the stimulated states, considering set tests for association (mtSet), interaction (iSet) and heterogeneity-GxC (iSet-het). For comparison we also applied a single-variant tests using the same variants, testing for association (mtLMM) and stimulus interaction (mtLMM-int). For single-variant tests, we estimated gene-level significance using the P value of the lead cis variant (adjusted within cis regions using eigenMT, [30]). Empirical P values for iSet and iSet-het were estimated from 30 parametric bootstraps per-gene and stimulus, combining all null LLRs across probes (resulting in 380,310 null LLRs per stimulus overall). Empirical P values for mtSet were obtained using the same permutation procedure as in [12]. Results from all methods were adjusted for multiple testing across probes using the Benjamini Hochberg procedure applied to each stimulus context separately. Reported results correspond to significant effects at genome-wide FDR < 5% (, ).

Best linear unbiased predictor from single-context set test

To illustrate the properties of the heterogeneity-GxC QTLs detected by iSet-het, we additionally considered univariate set tests in the same cis regions, however independently modeling each cellular context. At FDR<5% this analysis revealed 4,187, 4,786, 4,240 and 4,620 probes with an eQTL respectively in the naive, IFN-gamma, LPS2h and LPS24h states (). To estimate the cis-genetic contribution to gene-expression in each context we calculated the Best Linear Unbiased Predictor (BLUP) from the model as , where μ is the estimated mean, is the estimated variance explained by cis variants, is the cis realized relatedness matrix, −1 is the inverse of the total estimated covariance and is the gene-expression vector in the corresponding context.

Single-variant forward selection LMM

We used a single-variant forward selection LMM [33] to characterize eQTLs with significant heterogeneity-GxC effects. The model was fit considering up to three steps for gene and context, iteratively accounting for lead variant as additional fixed effect covariates when significant (FDR<5%). For each cellular context, region-based P values were adjusted for multiple testing across probes using the Benjamini Hochberg procedure for each of the three steps (only across probes that were tested at that specific step). This analysis yielded 15,756, 2,690 and 457 instances (across all genes and contexts) with one, two or three associations respectively ().

Results from step-wise analysis were used to annotate probes with significant heterogeneity-GxC. We denoted the 1,449 probes that have significant marginal associations in both contexts and independent lead variants (r2 < 0.20) as a shift in lead variants between the two contexts. Probes with a shared lead eQTL (significant in both contexts, lead variants r2 > 0.80) were annotated using secondary effects. Among the 4,186 probes with shared main effects, this analysis revealed context-specific secondary QTLs were identified for 999 genes. Context-specific secondary effects were defined when either i) the secondary effect was significant in only one of the two contexts or ii) the secondary effects lead variants were in low LD (r2 < 0.20) ().

Annotation of opposite-effect eQTLs

We classified the 4,186 eQTLs with shared lead eQTL into directionally consistent and opposite-effect eQTLs. Briefly, opposite effects were defined by three criteria, i) marginal significance in both contexts, ii) LD between contexts (r2>0.8) and iii) negative correlation of genetic effects. These criteria resulted in 67 opposite-direction QTLs. Directionally consistent eQTLs correspond to criteria i) and ii) but positive correlated genetic effects, resulting in 4,119 co-located QTLs. Statistical significance of the enrichment for significant heterogeneity-GxC effects in opposite-direction eQTLs rather than same-direction eQTLs was assessed using a one-sided Fisher’s exact test ().

iSet for analysis of stratified cohorts

Simulations for analysis of stratified individuals

To study performance of iSet when considering interaction analyses in stratified cohorts, we considered simulation experiments analogous to those for fully observed designs. We generated a synthetic cohort of 2,000 Europeans where each individual was phenotyped in only in one of two contexts. For each individual, the phenotyped context was independently selected using a draw from a Bernoulli distribution (symmetric, 50% success rate). Statistical calibration and power simulations were performed analogously to the approach used for fully observed designs. Population structure was accounted for using the first ten principal components of the realized relatedness matrix as fixed effect covariates. We did not consider tests for heterogeneity-GxC, as differential tagging of causal variants could potentially result in spurious heterogeneity-GxC signals, and hence additional controls would be required. However, in principle the test applies to stratified populations.

Comparison to alternative methods

We compared iSet to the single-variant interaction tests as in [10] (mtLMM-int) and the gene-environment set association test (GESAT) [13]. The latter approach is representative for a family of closely related set tests that can only be applied to test for interaction effects in stratified populations (See ). As an additional comparison, we extended the single-variant interaction test in [10] for stratified cohorts. To the best of our knowledge there are currently no implementations of mtLMM-int that can be applied to such designs. The models are available within the LIMIX package [48] (for full details see ). GESAT was run using the function GESAT of the package iSKAT version 1.2. Both iSet and GESAT were applied on identically processed standardized variants.

Genotype-sex interaction tests in lipid traits

We performed a genotype-sex interaction analysis of four blood lipid phenotypes (C-reactive protein (CRP), triglycerides (TRIGL), LDL and HDL cholesterol levels) measured in 5,256 unrelated individuals from the NFBC1966 cohort [20] (phs000276.v1.p1). Following [11, 12], we regressed out major covariates, following a quantile-normalization of each trait individually. In order to correct for population structure, we considered the first ten principal components of the realized relatedness matrix as fixed effect covariates.

We applied mtSet and iSet to 318,653 genome-wide variants with an allele frequency of at least 1% using a sliding-window approach (100kb regions, 50kb step size; resulting in 52,819 windows overall; (). For comparison we considered the single-variant interaction test [10], GESAT [13] and stSet [8], a univariate set test without stratification by sex. For each window we considered 100 permutations for mtSet and stSet and 100 parametric bootstraps for iSet and combined the obtained null LLRs across windows and traits (for a total of 21,127,600 null LLRs per test) to obtain empirical P values. Significance of the considered statistical tests was assessed at FWER = 10%. Summary results from all considered methods are reported in .

Imputation of NFBC1966 genotypes

Genotype data from NFBC1966.phs000276.v1.p1 were imputed using the 1000 Genomes Project phase 3 reference panel as described in the following. After aligning the dataset to the reference panel, we ran shapeit v2.r727 [50] with recommended parameters on each chromosome to produce haplotype estimates. We used impute2 v2.3.2 [51] with recommended parameters to impute untyped genotypes. Imputation was performed on chunks of approximately 5Mb. We merged region with less than 200 SNPs and avoided considering regions that span the centromere.

Supplementary methods.

Derivation and implementation details of the gene-context interaction set tests.

(PDF)

Click here for additional data file.

Type-1 error estimates on simulated data.

Shown are empirical type-1 error estimates for increasingly stringent significance level thresholds. Persistent genetic effects (No GxC effects) were simulated using the standard simulation parameters (a, parameters in but without simulating rescaling) and the same setting while considering a single causal variant (b), outlying samples (c) and epistatic interactions between randomly selected pairs of variants (d). These results show that iSet and iSet-het yield P values with controlled type-I error rates under different types of model misspecification. The corresponding QQ plots are shown in (main) and .

(PDF)

Click here for additional data file.

Simulation settings.

Simulated phenotype data were generated as sum of effects from variants in a designated causal region, effects from a relatedness component/population structure, effects from K = 10 unmeasured hidden confounding factors, and iid observation noise (). We fixed the variance explained by the region (v = 2%), and the fraction of shared background signal (α = 0.6), the fraction of residual variance that is explained by the hidden factors (β = 0.5). We considered variable numbers of causal variants and altered the extent of rescaling-GxC and heterogeneity-GxC. (a) The contributions to phenotypic variance of all simulated effects. (b,c) Parameter values considered for simulations of rescaling-GxC and heterogeneity-GxC effects, respectively. Individual parameters in (b,c) were varied one by one, while keeping others at their default values (bold face). When simulating general-GxC, and when varying the number of causal variants, the region effect correlation r was constrained to 0.2 < r < 0.8. See for full details.

(PDF)

Click here for additional data file.

Relationship between the proportionality factor of the effect sizes used in simulations, fold change and relative direction of genetic effects across contexts.

Shown is the relationship between the proportionality factor of the effect sizes in the two simulated contexts (η, x-axis in ), the corresponding fold change of the effect sizes, the fold change of the variance explained by the region and the relative directionality of the genetic signals. For −1 < η < 1, the absolute fold change increases for small absolute values of η and tends to infinity in the limit η → 0. The range η < 0 corresponds to genetic effect with opposite effects. The setting η = 0 corresponds to a local genetic effect that is specific to one context (no effect in the second context). Finally, η = 1 corresponds to no interaction (the fold change is zero and the direction of genetic effects is the same).

(PDF)

Click here for additional data file.

Tabular summary of results from the monocyte gene expression analyses.

The result table is provided as separate supplementary information file.

(XLSX)

Click here for additional data file.

Enrichment analysis of heterogeneity-eQTLs in opposite direction QTLs.

Breakdown of probe/stimulus pairs with shared lead variants, stratified by concordance of the effect direction (opposite-direction versus same-direction eQTLs) and significance of the heterogeneity-GxC test (heter vs No heter) in naive/IFN (a), naive/LPS2 (b) and naive/LPS24 (c).

(PDF)

Click here for additional data file.

Tabular summary of the gene-by-sex interaction analysis in human blood lipid traits from NFBC1966 cohort.

The result table is provided as separate supplementary information file.

(XLSX)

Click here for additional data file.

Computational complexity of iSet.

Shown is the computational complexity of iSet for alternative designs (either complete or stratified data designs), strategies to adjust for confounding (either principal components or a random effect) and analysis settings (number of variants R lower than the number of individuals N or not). Here, N denotes the number of individuals, R is the number of variants in the region, NPC is the number of principal components and t is the number of function evaluations of the optimizer. Operations that do not include t as factor refer to the cost of upfront computations that do not need to be evaluated in every step of the optimization. Operations that entail substantial computational burden (cubic operation in N) are highlighted in red.

(PDF)

Click here for additional data file.

Computational cost of iSet for alternative designs and cohort sizes.

Shown is the average CPU time (in seconds) for one in interaction set test using a 30kb region for alternative designs and cohort sizes. By default iSet uses principal components to adjust for confounding (PC). Alternatively, iSet can be combined with an additional random-effect (RE) to adjust for structure such as relatedness (see Methods). Reported CPU times are empirical averages to perform the tests for association (mtSet), interaction (iSet) and heterogeneity-GxC (iSet-het) and include the cost of 30 parametric boostraps for each test to estimate P values (Methods). CPU times are averaged across 100 genomic regions (with a size of 30 kb, ). Note that the reported costs do not include the up-front cost of the principal components and the up-front eigenvalue decomposition of the global relatedness matrix for the random-effect method. Runtime estimates were obtained using a single core of an Intel Xeon CPU E5-2670 2.60 GHz processor.

(PDF)

Click here for additional data file.

Alternative study designs supported by iSet.

iSet supports efficient interaction set tests both for complete designs (where each individual is phenotyped in all analyzed contexts, (a), and stratified cohorts (where each individual is phenotyped in only one of the analyzed contexts, (b).

(PDF)

Click here for additional data file.

Statistical calibration of iSet-het when only rescaling-GxC effects are simulated.

Shown is the QQ plot for the P values obtained from the heterogeneity-GxC interaction test (iSet-het) when simulating rescaling-GxC (without heterogeneity-GxC, Methods) for three different scenarios: (a) positive proportionality factor of effect sizes (0.3) and 1 causal variant, (b) positive proportionality factor of effect sizes (0.3) and 4 causal variants, (c) negative proportionality factor of effect sizes (-0.5) and 4 causal variants.

(PDF)

Click here for additional data file.

Assessment of calibration under different types of model mismatch.

Shown are the QQ plots for the P values obtained from the interaction set test (iSet) and the heterogeneity-GxC set test (iSet-het) when simulating different types of model mismatch. (a, b) Simulated phenotypes that violate the assumption of the infinitesimal model iSet is based on (which assumes an additive linear model of all variants in the set): i) a single causal variant in the region (a) and ii) epistatic interactions (GxG, b, see Methods). (c, d) Violations of Gaussian distributed residuals, by simulating outlying samples (Methods). Notably, while iSet and iSet-het using unnormalised phenotypes yield inflated P values (c), both methods are well calibrated after quantile normalization to a unit variance Gaussian distribution (d, Methods). All the results from the real data applications are based on quantile normalized phenotypes (Methods).

(PDF)

Click here for additional data file.

Distribution of the number of variants, the number of effective tests estimated by eigenMT and the average squared correlation within the testing regions in the different datasets.

From left to right: distribution of the number of variants across the analyzed regions; distribution of the number of effective tests as estimated by eigenMT; distribution of the average pairwise squared Pearson correlation (r2) across all variants in each region; scatter plot of the number of effective tests versus the number of variants. Shown in color is the within-region average correlation (r2) across all pairs of variants. From top to bottom: (a) 10,000 30kb regions from the simulated data based on 1000 Genomes individuals, (b) 100kb cis regions (centered on the TSS) considered in the cis stimulus eQTL analysis (288 individuals) and (c) 100kb regions considered in genotype-sex interaction analysis in the NFBC1966 cohort (5,402 individuals).

(PDF)

Click here for additional data file.

Simulation results for synthetic genotypes without LD.

Shown are results analogous to those presented in , however, considering synthetic genotypes without LD.

(PDF)

Click here for additional data file.

Comparison of iSet-het and single-variant strategies for discriminating rescaling from heterogeneity-GxC.

Receiver operating curve (a) and precision-recall curve (b) for alternative approaches to classify heterogeneity-GxC. Considered was the iSet test to score the extent of heterogeneity (iSet-het) and a baseline approach based on single-trait single-variant LMMs. Briefly, for the baseline model the considered score is defined as 1 − r2 (where r is the Pearson correlation coefficient between lead variants identified in each context) for regions with significant associations in both contexts (P-value thresholds 0.5, 0.01, 1e-3, 1e-4). Regions that were not marginally significant in either of the two contexts were assigned a score of zero (Methods).

(PDF)

Click here for additional data file.

Power of iSet and iSet-het when simulating heterogeneity-GxC effects and increasing numbers of causal variants.

Shown is the power of iSet, iSet-het and a single-variant interaction test (mtLMM-int) for detecting GxC effects when simulating heterogeneity-GxC for increasing numbers of causal variants (0 < r < 0.8, where r is the correlation of the simulated genetic effects in the two contexts).

(PDF)

Click here for additional data file.

Assessment of genetic variance estimates from iSet using different covariance models.

Shown are the estimates of the genetic variance explained by the set component across all simulated settings when considering different covariance models. The variance component estimates are from a model that either considers a full-rank covariance (full, general case), a rank-one covariance (rank1, only rescaling-GxC) and a block covariance matrix (block, which models only persistent genetic effects). Both designs with fully observed cohorts (complete—1,000 individuals and 2 contexts for a total of 2,000 trait measurements) and stratified samples (stratified—2,000 individuals and 2 contexts for a total of 2,000 trait measurements) are considered. While the full-rank iSet model yields calibrated variance components, other methods yield biased estimates in some settings. In particular, we considered scenarios with either rescaling-GxC effects (where we varied the number of causal SNPs and the proportionality factor of the effect sizes across the two contexts) or heterogeneity-GxC (where we vary the number of SNPs). For each simulated scenario, we considered 1,000 simulated regions and altered the variance explained by the region (we consider the values 2%, 5% and 10%). Grey horizontal lines denote the true simulated local genetic variance.

(PDF)

Click here for additional data file.

Power of iSet and iSet-het in sliding window experiments with different sizes of the testing regions.

(a,b) Power of iSet when simulating causal regions of different sizes (x-axis) and for sliding window analyses with increasing window sizes (y-axis) (see Methods). Considered were both settings with pure rescaling (a) and settings with heterogeneity-GxC effects (b). For comparison, we also considered the power of mtLMM-int (adjusted for multiple testing within the same windows using eigenMT, Methods). (c) Power of iSet in the same settings as considered in (b) when simulating heterogeneity-GxC effects and for the iSet-het test. Note that mtLMM-int is not specific to detecting heterogeneity-GxC; hence the reported power was set to zero.

(PDF)

Click here for additional data file.

Number of positives for single-variant methods and set tests as a function of the false discovery rate in the monocyte stimulus QTL data.

Shown is the number of probe/stimulus pairs with significant effects, detected by alternative single-variant and set-based tests, varying the genome-wide false discovery rate (FDR) threshold. Considered were the single-variant tests for associations (mtLMM-any) and interactions (mtLMM-int), as well as set tests for associations (mtSet), interactions (iSet) and heterogeneity-GxC effects (iSet-het). Venn diagrams on top show the overlap of significant probe/stimulus pairs for the three interaction tests at selected FDR thresholds (FDR<1%,2%,5%,10%).

(PDF)

Click here for additional data file.

Comparison of single-variant methods and set tests in the monocyte stimulus eQTL data.

Shown are the scatter plots of the −log10P values from single-variant LMMs and set tests for association tests (mtLMM-any vs mtSet) and interaction tests (mtLMM-int vs iSet) for different stimulus contexts (IFN/naive, LPS-2h/naive, LPS-24h/naive). P values for single-variant models correspond to the minimum P-value across variants in the considered testing region, adjusted for the effective number of tests (estimated using eigenMT, Methods). Venn diagrams on the top of individual panels show the overlap of probes with significant associations or interactions identified using alternative methods (5% FDR).

(PDF)

Click here for additional data file.

Results from single-trait set tests applied to individual cellular contexts.

(a) Scatter plot of the maximum Q values across pairs of jointly analyzed contexts (independent analysis), comparing set tests and single-variant tests for different stimulus pairs. Heterogeneity-GxC cases without clear single-variant interpretation (highlighted in red) tend to be associated with increased significance when using set tests, suggesting that differences in power may impair the single-variant annotation of these eQTLs. (b) Histogram of the squared Pearson correlation (r2) of cis genetic effects predicted using the best linear unbiased predictor, BLUP (Methods) across pairs of contexts. Shown is the distributions either for all probes with significant effects in the naive context and probes with significant heterogeneity-GxC. As expected, genes with heterogeneity-GxC tend to be associated with lower correlation of genetic effects across contexts. (c) Scatter plot of the correlation coefficient of genetic effects (as in b) versus the maximum q-value across the naive and the stimulated state. Heterogeneity eQTLs (highlighted in red) tend to have stronger marginal associations and lowly-correlated cis genetic signals.

(PDF)

Click here for additional data file.

Examples of opposite-effect eQTLs with significant heterogeneity-GxC effects.

Shown is the z-score statistics for cis variants for OAS (a), LMNA (b) and PTK2B (c) across the contexts showing opposite effects. While the three examples are identified as opposite effects when using single-variant methods (Fairfax et al., Science, 2014), iSet identified significant heterogeneity-GxC, suggesting changes in the configuration of causal variants. Lead variants in individual contexts are annotated using triangles and are in high LD (r2 > 0.8).

(PDF)

Click here for additional data file.

Calibration and power simulations for different interaction methods for analysis of stratified cohorts.

(a) QQ plot for the P values obtained when applying iSet and GESAT to synthetic datasets where only persistent genetic effects (No GxC) were simulated. (b) Power comparison of alternative interaction tests and variance decomposition results from iSet when simulating rescaling-GxC effects, for different factors of proportionality of the variant effect sizes across contexts. Considered were iSet, GESAT and a single-variant interaction test (mtLMM-int). (c) Analogous results as in b when simulating general GxC effects (both rescaling-GxC and heterogeneity-GxC). Results are stratified by the correlation between the simulated genetic effects between the two contexts.

(PDF)

Click here for additional data file.

Manhattan plots from alternative methods applied for the genome-wide analysis of human lipid levels in the NFBC1966 cohort.

Shown are Manhattan plots for C-reactive protein (crp3dec, a), LDL cholesterol (FS_KOL_L, b), HDL cholesterol (FS_KOL_H, c), and triglycerides (FS_TRIG, d) obtained from univariate set tests ignoring sex-specific differences (stSet), an association test that accounts for differences in genetic effect across strata (mtSet), iSet, GESAT and single-variant interaction test (mtLMM-int). Red arrows indicate the interaction effects that are discussed in the main text. Blue arrows indicated associations that can only be detected when modeling differences in effect sizes across strata (mtSet vs stSet).

(PDF)

Click here for additional data file.

| QQ plots when applying alternative methods to lipid levels in NFBC1966.

Shown are the QQ plots for C-reactive protein (crp3dec), LDL cholesterol (FS_KOL_L), HDL cholesterol (FS_KOL_H), and triglycerides (FS_TRIG) obtained using a univariate association set tests ignoring sex (stSet, a), an association test modeling sex-specific genetic effects (mtSet, b), iSet (c), GESAT (d) and single-variant interaction test (mtLMM-int).

(PDF)

Click here for additional data file.

Manhattan plot in the interaction locus for C-reactive protein using single-variant interaction tests on imputed variants.

Shown is the Manhattan plot for C-reactive protein using single-variant interaction tests applied to common variants (MAF>0.5%) on imputed data. Vertical grey lines indicate the 100kb region with significant genotype-sex interaction (FWER = 10%) detected using iSet. Non-imputed typed variants are highlighted in red, showing that for this locus imputation strategies did not increase the power of single-variant methods.

(PDF)

Click here for additional data file.

The interaction for C-Reactive protein on chromosome 1 is a male-specific effect.

(a) Local Manhattan plot (1Mb around significant region) for single-variant association tests, either considering males (black) or females (pink). For comparison, shown are also the P values from the iSet (red), mtSet (green) and stSet (grey). (b) C-Reactive protein level stratified by different alleles of rs4660378 (lead SNP identified in the analysis using male individuals only). Rho in the caption correspond to the Spearman rank correlation coefficient and the corresponding P value, both for male and female individuals.

(PDF)

Click here for additional data file.