Kernel-based gene–environment interaction tests for rare variants with multiple quantitative phenotypes

Previous studies have suggested that gene–environment interactions (GEIs) between a common variant and an environmental factor can influence multiple correlated phenotypes simultaneously, that is, GEI pleiotropy, and that analyzing multiple phenotypes jointly is more powerful than analyzing phenotypes separately by using single-phenotype GEI tests. Methods to test the GEI for rare variants with multiple phenotypes are, however, lacking. In our work, we model the correlation among the GEI effects of a variant on multiple quantitative phenotypes through four kernels and propose four multiphenotype GEI tests for rare variants, which are a test with a homogeneous kernel (Hom-GEI), a test with a heterogeneous kernel (Het-GEI), a test with a projection phenotype kernel (PPK-GEI) and a test with a linear phenotype kernel (LPK-GEI). Through numerical simulations, we show that correlation among phenotypes can enhance the statistical power except for LPK-GEI, which simply combines statistics from single-phenotype GEI tests and ignores the phenotypic correlations. Among almost all considered scenarios, Het-GEI and PPK-GEI are more powerful than Hom-GEI and LPK-GEI. We apply Het-GEI and PPK-GEI in the genome-wide GEI analysis of systolic blood pressure (SBP) and diastolic blood pressure (DBP) in the UK Biobank. We analyze 18,101 genes and find that LEUTX is associated with SBP and DBP (p = 2.20×10−6) through its interaction with hemoglobin. The single-phenotype GEI test and our multiphenotype GEI tests Het-GEI and PPK-GEI are also used to evaluate the gene–hemoglobin interactions for 22 genes that were previously reported to be associated with SBP or DBP in a meta-analysis of genetic main effects. MYO1C shows nominal significance (p < 0.05) by the Het-GEI test. NOS3 shows nominal significance in DBP and MYO1C in both SBP and DBP by the single-phenotype GEI test.


Introduction
Genome-wide association studies (GWASs) have identified numerous common variants associated with common diseases or phenotypes [1]. Nevertheless, a small portion of the heritabilities can be explained by the discovered common variants [2,3]. Sequencing studies showed that some of the "missing heritability" was attributable to rare variants [4,5]. Complex  are usually influenced by genetic factors, environmental factors and the interplay between them. Wang et al. showed that the interactions between SMC5 variants and alcohol consumption are associated with fasting plasma lipid levels [6]. Yang et al. demonstrated that the interactions between PDE3B variants and smoking are associated with pulmonary function [7]. Johansson et al. revealed that the interactions between NFE2L2 variants and second-hand smoke are associated with pediatric asthma risk [8]. For a long time, gene-environment interactions (GEIs) have been expected to explain some of the "missing heritability" and shed light on the genetic etiology of complex diseases [9]. Existing studies suggest that the interaction between a common variant and an environmental factor may be associated with multiple correlated phenotypes, which is called GEI pleiotropy [10]. Kilpeläinen et al. identified four loci in or near CLASP1, LHX1, SNTA1 and CNTNAP2 that are associated with three blood lipid levels: low density lipoprotein, high density lipoprotein and triglycerides through their interactions with physical activity [11]. Novel gene-sleep interactions were also identified for known lipid loci, including LPL and PCSK9 [12]. To date, all the reported GEI pleiotropies are with common variants. From a methodological perspective, Majumdar et al. showed that statistical power to detect GEI effects can be improved by analyzing multiple phenotypes jointly [10]. However, multiphenotype methods for testing GEIs with rare variants are lacking.
To the best of our knowledge, there is only one method currently available for testing GEIs with rare variants and multiple phenotypes [13]. The method consists of three steps: remove correlation among multiple phenotypes by using principal component analysis or other linear transformations; obtain p value for each transformed phenotype by testing the effects of an optimally weighted combination of GEIs for rare variants (TOW-GE) [14]; employ Fisher's combination test (FCT) to combine the p values of multiple phenotypes. We denote the method as TOWGE-FCT in this paper. It can be expected that the degree of freedom of TOW-GE-FCT would become larger with the increasing number of phenotypes, which might limit statistical power of the test.
In this work, we model the correlations among the GEI effects of a variant on multiple phenotypes by assuming four different kernel matrices, similar to those for multiphenotype tests of genetic main effects [15]. We extend the single-phenotype GEI test [16] and propose four multiphenotype GEI tests for rare variants, which are the test with homogeneous kernel (Hom-GEI), the test with heterogeneous kernel (Het-GEI), the test with projection phenotype kernel (PPK-GEI) and the test with linear phenotype kernel (LPK-GEI). We conduct simulation studies to examine the empirical distributions of the four test statistics under the null hypothesis and compare their statistical power under different scenarios. In the analysis of systolic blood pressure (SBP) and diastolic blood pressure (DBP) in the UK Biobank, we chose hemoglobin (Hb) as the environmental variable, which is known to be associated with both SBP and DBP [17,18]. With the whole-exome sequencing data in 200,643 samples, we applied Het-GEI and PPK-GEI in the genome-wide analyses of gene-Hb interactions. We also carried out single-phenotype and multiphenotype GEI tests to evaluate the gene-Hb interactions for 22 genes that were previously reported to be associated with SBP or DBP in a meta-analysis of main genetic effects [19].

Single-phenotype GEI test
Assume that n unrelated individuals are sequenced in a gene or region with m rare variants and K quantitative phenotypes are measured. For the k-th phenotype, y k = (y 1k , y 1k , � � �, y nk ) T denotes an n × 1 phenotype vector, and X = (X 1 , X 2 , � � �, X q+1 ) is an n × (q + 1) matrix comprised of intercept and covariate vectors with X t = (X 1t , X 2t , � � �, X nt ) T , t = 1, 2, . . ., q+1. The first vector X 1 represents the intercept vector with elements X i1 = 1 and i = 1, 2, � � �, n. The other q vectors are the covariate vectors. Let G = (G 1 , G 1 , � � �, G m ) be an n × m genotype matrix, in which G j = (G 1j , G 2j , � � �, G nj ) T , j = 1, 2, . . ., m, and G ij is the number of minor alleles. E = diag {E i } denotes an n × n diagonal matrix of environmental measurements, and E i is centralized and included in X as a covariate for adjusting the environmental effect. Following the single-phenotype GEI test for rare variants in rareGE [16], we consider the linear mixed model as follows: where α k = (α k1 , α k2 , � � �, α k(q+1) ) T is a (q+1)×1 vector of covariate effects for the k-th phenotype. W = diag{w j } is an m × m weight matrix for the m variants. The weight of the j-th variants is w j = Beta(MAF j , 1, 25) [20], where MAF j is the minor allele frequency (MAF) of the j-th variants.
In addition, β k = (β 1k , β 2k , � � �, β mk ) T is an m × 1 vector consisting of genetic main effects for the k-th phenotype, and γ k = (γ 1k , γ 2k , � � �, γ mk ) T is an m×1 vector of the interaction effects. Here, the main genetic effects β k are assumed to be fixed and the interaction effects γ k to be random, γ k * MVN(0, σ 2 I m ). In addition, ε k = (ε 1k , ε 2k , � � �, ε nk ) T denotes an n×1 error vector, and ε k � MVNð0; s 2 k I n Þ. The null hypothesis for testing the GEI interactions is H 0 : σ 2 = 0. The model under the null hypothesis is Here, α k , β k , and s 2 k can be estimated by linear regression, and the estimated mean and variance-covariance matrix of y k areμ GW). The score statistic for testing the GEI effects is which is mathematically equivalent to Here, S jk ¼ P n i¼1 E i G ij ðy ik Àm ik Þ=ŝ 2 k is the score statistic for the j-th variant. Under H 0 , Q k � P j l j w 2 1;j follows a mixture of chi-square distributions with 1 degree of freedom, and λ j are nonzero eigenvalues of the regional relationship matrix The p-value can be computed by using Kuonen's saddlepoint approximation method [16,21]. In the same spirit, we extend the single-phenotype GEI test for multiple phenotypes.

Kernel-based multiphenotype GEI tests
Denote y = (y 1 , y 2 , � � �, y K ) as the n × K matrix of K phenotypes and A = (α 1 , α 2 , � � �, α K ) as the (q+1) × K matrix of covariate effects. Let B = (β 1 , β 2 , � � �, β K ) be the m × K matrix of genetic main effects and Γ = (γ 1 , γ 2 , � � �, γ K ) be the m × K matrix of GEI effects. In addition, ε = (ε 1 , ε 2 , � � �, ε K ) is the n × K error matrix. In light of the correlation among phenotypes, we assume ε = (ε 1 , ε 2 , � � �, ε K ) * MVN(0, Σ), i = 1, 2, � � �, n. Then, the mixed model for multiple phenotypes can be formulated in a matrix form as follows: Stack columns of the phenotype matrix y into a vector vecðyÞ ¼ ðy T 1 ; y T 2 ; � � � ; y T K Þ T and columns of the error matrix ε into vecðεÞ ¼ ðε where � is the Kronecker product [22]. We rewrite model (6) in vector form as Assume vec(Γ)~MVN(0, σ 2 Σ P � I m ), where Σ P is a K × K kernel in the phenotype space and models the correlation among the GEI effects of a variant on multiple phenotypes. As a result, vec(y)~MVN(vec(μ), H), where μ = XA + GWB and H = σ 2 (Σ P � EGWWG T E) + Σ � I n . The null hypothesis for testing the GEI effects with multiple phenotypes is H 0 : σ 2 = 0, and the score statistic is whereμ andΣ are the estimated mean and variance-covariance matrix, respectively. The score statistic Q asymptotically follows a mixture of 1-freedom chi-square distributions X j l j w 2 1;j and λ j are nonzero eigenvalues of where P = I n −Z(Z T Z) −1 Z T and Z is the same as in the single-phenotype GEI test. The corresponding p-values can be computed via Kuonen's saddlepoint method [21]. As can be seen in (8), our proposed tests depend on the kernel matrix Σ P . Similar to [15], we use four types of kernel matrices to model the correlation among the GEI effects on multiple phenotypes.
Homogeneous kernel. Assume that the GEI effects of a variant on multiple different phenotypes are homogeneous, implying that γ j1 = γ j2 = � � � = γ jK . The kernel is constructed as vector. Σ Hom indicates the GEI effects of a variant on multiple phenotypes to be the same. Heterogeneous kernel. Assuming that the GEI effect sizes of a variant on multiple phenotypes are heterogeneous, the kernel is where I K is a K × K identity matrix. Here, Σ Het implies that the GEI effects of a variant on multiple phenotypes are independent.
Projection phenotype kernel. Assume that the correlation among the GEI effects of a variant on multiple phenotypes can be depicted by the correlation among the phenotypes. That is, whereΣ is the estimated variance-covariance matrix of the phenotypes.
Linear phenotype kernel. Assume that the GEI effects of a variant on multiple phenotypes equal the squared correlation among the phenotypes. That is, Similar to the proof in [22], the test score statistic (8) is which can be rewritten as Therefore, the LPK-GEI test simply combines statistics of single-phenotype GEI tests across multiple phenotypes.
Based on different choices of the kernel matrix Σ P , we propose four multiphenotype GEI tests, which are named Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI.

Numerical simulations
To evaluate the null distributions and statistical power of the four proposed tests, we carried out extensive simulation studies. Using the calibrated coalescent model implemented in COSI [23], we generated 10,000 haplotypes in a 200 kb genomic region. Parameters in the coalescent model were used to mimic the linkage disequilibrium pattern, local recombination rate and demographic history for the population of European ancestry. We randomly paired these haplotypes to form diploid genotype data of 10,000 individuals and randomly selected 5000 out of the 10,000 individuals. A subregion length of 3 kb was randomly selected from the 200 kb region to obtain the genotype data of the 5000 samples for each replicate, and 1000 replicates of genotype data were generated. Variants with MAF � 0.01 were considered to be rare and used for simulations.
To evaluate null distributions of our proposed tests, four phenotypes of the 5000 unrelated individuals under the null hypothesis were generated. For the sake of simplicity, phenotypes shared the same covariate sets and were generated as follows: Where y i1 , y i2 , y i3 , y i4 are the four phenotypes for the i-th individual (i = 1, 2, . . ., n). For the i-th individual, sex i is a binary covariate following a Bernoulli distribution with probability 0.5, namely, sex i~B ernoulli(0.5). Both age i and bmi i are continuous covariates: age i~N (50, 25), bmi i~N (50, 25). G ij (j = 1, 2, � � �, m) are the coded genotypes of the simulated causal variants for individual i. Here, we assumed a proportion of causal variants θ = 0.1, 0.2, 0.3. In addition, w j (j = 1, 2, � � �, m) is the weight of variant j; β j1 , β j2 , β j3 , and β j4 are the main genetic effects for phenotype 1, phenotype 2, phenotype 3 and phenotype 4, respectively, with β j1 = 0.1, β j2 = 0.2, β j3 = 0.1, and β j4 = 0.2; and l 1 , l 2 , l 3 , and l 4 are indicator variables, with l k = 1 when phenotype k is associated with genetic variants and l k = 0 otherwise. Since not all of the phenotypes may be associated with the rare variants [22,24], we considered scenarios under which pleiotropy exists or does not exist. Specifically, we assumed that the first t phenotypes were associated with the rare variants, namely, l 1 = � � � = l t = 1 and l t+1 = � � � = l 4-t = 0, t = 1,2,3,4. ε ik is a random error for the i-th individual and the k-th phenotype, To evaluate the statistical power of our proposed four multiphenotype GEI tests, we simulated the four correlated phenotypes for 5000 independent individuals under the alternative hypothesis. For each of the genotype datasets, one phenotype and covariates set was simulated according to the following model: where sex i , age i , bmi i , G ij , w j , β jk (k = 1,2,3,4), l k (k = 1,2,3,4) and ε i = (ε i1 , ε i2 , ε i3 , ε i4 ) T are the same as described in model (12). The body mass index (BMI) was centered and used as the environmental variate E i . Here, γ jk is the gene-BMI interaction effect of the j-th causal rare variant on the k-th phenotype, with γ jk~N (0, 0.05) 2 . Since the interaction effects of a variant on each phenotype were simulated independently, the gene-BMI interaction effects of a variant on multiple phenotypes are heterogeneous. In all simulations and the analyses of the simulated data, variant weights were the density function of beta distribution with degrees of freedom of 1 and 25 evaluated at the MAF of rare variants [20] as described in the single-phenotype GEI test. We considered the gene-BMI interaction to be significant if its p-value was less than 2.5×10 −6 , corresponding to a correction for multiple testing in genome-wide studies of 20,000 genes. Empirical power was the portion of significant results in 1000 replicates.
Null distributions. We examined null distributions of test statistics for Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI with causal rare variant proportion θ = 0.2 and all phenotypes associated with genetic variants. We first estimated means and residuals by performing phenotype-specific regression analyses. The variance-covariance matrix Σ was estimated by residuals from all phenotypes. Test statistics of the four tests were computed as in (8) with corresponding kernels for Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI. Using Kuonen's saddlepoint approximation method [21], the p-values of the four test statistics were computed. Finally, we compared distributions of empirical p-values with the expected uniform distribution between 0 and 1.
The quantile-quantile (Q-Q) plots of the four multiphenotype statistics under weak (ρ = 0.25), moderate (ρ = 0.5) and strong (ρ = 0.75) correlations among phenotypes are shown in From Figs 4-6, we can see that the four tests provide improved power with more phenotypes associated with the interactions, which suggests that our multiphenotype analyses can  Fig 6C. This demonstrates that Hom-GEI, Het-GEI and PPK-GEI can benefit from the increased correlations among phenotypes. However, since LPK-GEI directly combines statistics from singlephenotype GEI tests and ignores the phenotypic correlations, the increased correlation among phenotypes leads to a substantial power loss.
Among almost all of the considered scenarios in Figs 4-6, PPK-GEI has approximately the same or slightly larger power than Het-GEI, and the two tests outperform Hom-GEI and LPK-GEI. Hom-GEI shows the poorest power performance among all the proposed tests. This is because the phenotypes were simulated based on heterogeneous interaction effects, violating the assumption that Hom-GEI is based upon. Because the GEI effects of a variant on multiple phenotypes can hardly be homogeneous in reality, Hom-GEI may not be a good choice for real data analysis. Therefore, we choose Het-GEI and PPK-GEI for our genome-wide interaction analysis in UK Biobank.
With the proportion of causal rare variant θ = 0.1 and the among-phenotype correlation ρ = 0.5, we compared the power of our tests with the TOWGE-FCT under different numbers of phenotypes associated with the interactions, the results are shown in Fig 7. Because TOW-GE-FCT is a permutation based method, which is computationally very expensive, the power results were evaluated at the significance level of 0.05. As can be observed from Fig 7 that all the five tests can provide enhanced power as more phenotypes associated with the interactions. Het-GEI, PPK-GEI and LPK-GEI tests have higher power than TOWGE-FCT, however, Hom-GEI has lower power than TOWGE-FCT. We can also see that Het-GEI and PPK-GEI tests outperform the other tests, further indicating Het-GEI and PPK-GEI to be two powerful tests. For instance, when only the first two phenotypes are associated with the interactions, the power values for Hom-GEI, Het-GEI, PPK-GEI, LPK-GEI and TOWGE-FCT are 0.177, 0.507, 0.517, 0.403 and 0.304, respectively.

Gene-Hb interaction analysis of blood pressure phenotypes in UK Biobank
UK Biobank is a prospective study that recruited approximately 500,000 volunteers aged 40 and 69 years in the United Kingdom and collected extensive genetic and phenotypic data [25,26]. We used the whole-exome sequencing data released by UK Biobank with a total of 200,643 samples. Individuals who had withdrawn and one member in each pair with kinship larger than 0.25 measured via KING [27] were removed. We considered SBP and DBP as the blood pressure (BP) phenotypes and Hb as the environmental factor, which is known to be associated with both SBP and DBP [17,18]. Covariates included age, age 2 , sex, BMI and 20 principal components to adjust for population stratification. SBP and DBP averaged over multiple measurements at baseline were used. For individuals taking BP-lowering medications, 10 mm Hg and 5 mm Hg were added to the SBP and DBP, respectively [28,29]. Phenotypes and covariates located 5 standard deviations away from their respective means were defined as

PLOS ONE
outliers. Outliers or individuals with missing phenotypes or missing covariates were removed. As a result, 157,514 individuals, including 71,501 males (45.4% males) and 86,013 females (54.6% females), were included in our gene-Hb interaction analyses. BP phenotypes, age, BMI and Hb were standardized before the analysis.
For each of the 18,101 genes, we performed Het-GEI and PPK-GEI analysis of the gene-Hb interactions on SBP and DBP phenotypes. Manhattan plots of p-values from the two tests are presented in Fig 8, and QQ plots of the Het-GEI and PPK-GEI tests are shown in Fig 9. With a genome-wide significance level of 2.5×10 −6 , only LEUTX is significant according to the Het-GEI test (p-value = 2.2×10 −6 ), and its p-value according to the PPK-GEI test is 7.43×10 −6 . If we

PLOS ONE
consider a suggestive significance level at 1×10 −4 , twelve genes passed the threshold, whose details are presented in Table 1.
Recently, Surendran et al. reported 22 genes associated with SBP or DBP in a meta-analysis of 1.3 million samples from multiple cohorts, including UK Biobank, the Million Veterans Program and deCODE [19]. For the 22 genes, we looked up our genome-wide results for the possible interactions with Hb. For comparison, we also conducted single-phenotype GEI tests using the INT-FIX function from the rareGE R package [16] for the two BP phenotypes

PLOS ONE
separately. The number of rare variants involved in the analysis and p-values from the multiphenotype and single-phenotype GEI tests are provided in Table 2.
There are no significant results after correcting the multiple testing. At the nominal significance level of 0.05, only one gene, MYO1C, shows interactions with Hb for BP phenotypes by the PPK-GEI test (p-value = 0.038). With the single-phenotype GEI test, NOS3 has a p-value of 0.026 for DBP, and MYO1C has p-values of 0.018 and 0.011 for SBP and DBP, respectively.

Discussion
In this paper, we propose four statistical tests, Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI, to test GEI effects with rare variants for multiple correlated quantitative phenotypes. Through simulation studies, the statistical power of the tests was investigated in terms of the proportion of causal variants, the number of phenotypes associated with interactions and the correlation strength among phenotypes. Simulation results show that all tests demonstrate improved statistical power when the proportion of causal variants or the number of associated phenotypes increases. Hom-GEI, Het-GEI and PPK-GEI benefit from correlation among phenotypes; however, the LPK-GEI test suffers power loss, especially when correlation among phenotypes is strong. This is because LPK-GEI directly combines statistics from single-phenotype GEI tests and ignores phenotype dependence. In addition, among almost all of the considered scenarios, Het-GEI and PPK-GEI have almost the same power and outperform the other two tests. Hom-GEI shows the poorest power due to its unrealistic assumption. In summary, Het-GEI and PPK-GEI are two powerful tests for investigating GEI with multiple quantitative phenotypes.
We applied Het-GEI and PPK-GEI in the genome-wide analysis of SBP and DBP in order to detect possible gene-Hb interactions in UK Biobank. We analyzed 18,101 genes and identified LEUTX to be associated with BP phenotypes through its interaction with Hb via the Het-GEI test. At the suggestive significance level, twelve genes were identified to be associated with BP phenotypes through their interactions with Hb. LEUTX was previously reported to play a central role in embryo genome activation [34] whose role in BP regulation is unclear. Recent study of rare variants suggests that BP-associated variants are enriched in active chromatin regions of fetal tissue and potentially link fetal development to BP regulation in later life [19]. Thus, LEUTX might be associated with the BP phenotypes in a similar manner. In the analysis of 454,787 UK Biobank participants, genetic main effect of LEUTX was not associated with BP phenotypes at the nominal significance level of 0.05 [35]. However, our tests identified LEUTX to be interacted with Hb in BP phenotypes at the genome-wide significance level of 2.5×10 −6 in the sample of 157,514 UK Biobank participants. We also conducted a single-phenotype GEI test and multiphenotype GEI tests to evaluate the gene-Hb interactions for 22 genes that were previously reported to be associated with SBP or DBP in a meta-analysis of genetic main effects. MYO1C shows nominal significance by the Het-GEI test. NOS3 shows nominal significance in DBP and MYO1C in SBP and DBP by the single-phenotype GEI test. Our proposed multiphenotype GEI tests are an extension of the single-phenotype GEI test in rareGE [16]. The tests retain the desirable properties of the single-phenotype GEI test, which allows for adjusting covariates and is powerful when the GEI effects of variants act in different directions on phenotypes. Our proposed multiphenotype GEI tests, except for Hom-GEI, are more powerful than the existing multiphenotype GEI test TOWGE-FCT. Moreover, our methods are computationally less expensive. This is because that TOWGE-FCT employs permutations to evaluate p value for each transferred phenotype, however, our tests are based on asymptotic distributions and p values can be computed analytically.
Our proposed multiphenotype GEI tests have the following limitations. First, our tests have specific assumptions on the correlation structure for the GEI effects among multiple phenotypes, violating the assumption would lead to a substantial loss of power. Second, our tests assume individuals to be unrelated. Familiar correlation is not considered, and related samples must be excluded, which may substantially reduce the sample size and result in a loss of power. Third, our proposed GEI tests are for rare variants only, while both common and rare variants may contribute to complex diseases [36][37][38]. In future work, we plan to consider a unified test to address these issues. Finally, while we identified LEUTX to interact with Hb in BP phenotypes, the result lacks independent validation. Thus, it should be considered as being preliminary and further experiments are necessary.

Conclusion
In this paper, we modeled the correlation among the GEI effects of a variant on multiple phenotypes by using four kernels. Based on these kernels, we proposed four multiphenotype GEI tests for rare variants. The four tests retain the desirable properties of the single-phenotype GEI test and provide enhanced statistical power by analyzing multiple phenotypes simultaneously. We applied Het-GEI and PPK-GEI to test gene-Hb interactions for 18,101 genes in SBP and DBP in UK Biobank. LEUTX was associated with BP phenotypes through the interaction with Hb via the Het-GEI test. At the suggestive significance level, twelve genes were reported. Our proposed tests can be readily used to test GEIs in a variety of correlated phenotypes and hopefully contribute to the genetic studies of complex diseases.