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 diseases 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 TOWGE-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].

Methods

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, = (y1, y1, ⋯, y)T denotes an n × 1 phenotype vector, and = (1, 2, ⋯, ) is an n × (q + 1) matrix comprised of intercept and covariate vectors with = (X1, X2, ⋯, X)T, t = 1, 2, …, q+1. The first vector 1 represents the intercept vector with elements X = 1 and i = 1, 2, ⋯, n. The other q vectors are the covariate vectors. Let = (1, 1, ⋯, ) be an n × m genotype matrix, in which = (G1, G2, ⋯, G)T, j = 1, 2, …, m, and G is the number of minor alleles. = diag{E} denotes an n × n diagonal matrix of environmental measurements, and E is centralized and included in 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 = (α, α, ⋯, α)T is a (q+1)×1 vector of covariate effects for the k-th phenotype. = diag{w} is an m × m weight matrix for the m variants. The weight of the j-th variants is w = Beta(MAF, 1, 25) [20], where MAF is the minor allele frequency (MAF) of the j-th variants. In addition, = (β1, β2, ⋯, β)T is an m × 1 vector consisting of genetic main effects for the k-th phenotype, and = (γ1, γ2, ⋯, γ)T is an m×1 vector of the interaction effects. Here, the main genetic effects are assumed to be fixed and the interaction effects to be random, ~ MVN(0, σ2). In addition, = (ε1, ε2, ⋯, ε)T denotes an n×1 error vector, and . The null hypothesis for testing the GEI interactions is H0: σ2 = 0. The model under the null hypothesis is

Here, , , and can be estimated by linear regression, and the estimated mean and variance-covariance matrix of are where , with = (, ). The score statistic for testing the GEI effects is which is mathematically equivalent to Here, is the score statistic for the j-th variant.

Under H0, follows a mixture of chi-square distributions with 1 degree of freedom, and λ 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 = (1, 2, ⋯, ) as the n × K matrix of K phenotypes and = (1, 2, ⋯, ) as the (q+1) × K matrix of covariate effects. Let = (1, 2, ⋯, ) be the m × K matrix of genetic main effects and = (1, 2, ⋯, ) be the m × K matrix of GEI effects. In addition, = (1, 2, ⋯, ) is the n × K error matrix. In light of the correlation among phenotypes, we assume = (ε1, ε2, ⋯, ε) ~ 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 into a vector and columns of the error matrix into . We have vec()~MVN(0, ⊗ ), where ⊗ is the Kronecker product [22]. We rewrite model (6) in vector form as

Assume vec()~MVN(0, σ2P ⊗ ), 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()~MVN(vec(), ), where = + and = σ2(P ⊗ ) + ⊗ . The null hypothesis for testing the GEI effects with multiple phenotypes is H0: σ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 and λ are nonzero eigenvalues of where = −(T)−1T and 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 γ = γ = ⋯ = γ. The kernel is constructed as where 1 = (1, 1, ⋯, 1)T is a K × 1 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 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.

Results

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, y, y, y are the four phenotypes for the i-th individual (i = 1, 2, …, n). For the i-th individual, sex is a binary covariate following a Bernoulli distribution with probability 0.5, namely, sex ~ Bernoulli(0.5). Both age and bmi are continuous covariates: age ~ N(50, 25), bmi ~ N(50, 25). G (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 = 1, 2, ⋯, m) is the weight of variant j; β, β, β, and β are the main genetic effects for phenotype 1, phenotype 2, phenotype 3 and phenotype 4, respectively, with β = 0.1, β = 0.2, β = 0.1, and β = 0.2; and l1, l2, l3, and l4 are indicator variables, with l = 1 when phenotype k is associated with genetic variants and l = 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, l1 = ⋯ = l = 1 and l = ⋯ = l4- = 0, t = 1,2,3,4. ε is a random error for the i-th individual and the k-th phenotype, = (ε, ε, ε, ε)T ~ MVN(0, ), where and ρ represents the correlation among different phenotypes. Three levels of correlation strength were considered: weak correlation with ρ = 0.25, moderate correlation with ρ = 0.5 and strong correlation with ρ = 0.75. For each simulation setup, 20 replicates of phenotypes and covariates were simulated based on one genotype dataset; thus, a total of 20,000 replicates of phenotypes and covariates were simulated.

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, age, bmi, G, w, β (k = 1,2,3,4), l (k = 1,2,3,4) and = (ε, ε, ε, ε)T are the same as described in model (12). The body mass index (BMI) was centered and used as the environmental variate E. Here, γ is the gene-BMI interaction effect of the j-th causal rare variant on the k-th phenotype, with γ ~ 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 Figs 1–3, respectively. The empirical distributions of the four test statistics are aligned with their theoretical distributions, as expected.

Fig 1

Q–Q plots of the test statistics under the null hypothesis with weak among-phenotype correlation ρ = 0.25.

The horizontal and vertical axes represent the negative log10 of the expected p-values and the negative log10 of the observed p-values, respectively. A: Hom-GEI; B: Het-GEI; C: PPK-GEI; D: LPK-GEI.

Fig 3

Q–Q plots of the null distributions under the null hypothesis with strong among-phenotype correlation ρ = 0.75.

The horizontal and vertical axes represent the negative log10 of the expected p-values and the negative log10 of the observed p-values, respectively. A: Hom-GEI; B: Het-GEI; C: PPK-GEI; D: LPK-GEI.

Q–Q plots of the test statistics under the null hypothesis with moderate among-phenotype correlation ρ = 0.5.

Statistical power

The statistical power of Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI under weak, moderate and strong correlations among phenotypes are shown in Figs 4–6, respectively. In each figure, the power with three different proportions of causal variants and different numbers of associated phenotypes are presented. All four tests show improved power as the proportion of causal variants increases. This is because an increased proportion of causal variants leads to larger interaction effects under the test. Taking Fig 5D as an example, for the causal rare variant proportion θ = 0.1, Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI have powers of 0.067, 0.424, 0.431 and 0.199, respectively. For θ = 0.2, the corresponding powers are 0.145, 0.72, 0.725 and 0.428, respectively. For θ = 0.3, the power increases further to 0.278, 0.867, 0.869 and 0.608.

Fig 4

Statistical power of Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI under weak among-phenotype correlation ρ = 0.25.

The horizontal and vertical axes represent the proportion of causal rare variants and the statistical power, respectively. A: Power when only the first phenotype is associated with interactions; B: Power when the first two phenotypes are associated with interactions; C: Power when the first three phenotypes are associated with interactions; D: Power when all four phenotypes are associated with interactions.

Fig 6

Statistical power of Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI under strong among-phenotype correlation ρ = 0.75.

The horizontal and vertical axes represent the proportion of causal rare variants and the statistical power, respectively. A: Power when only the first phenotype is associated with the interactions; B: Power when the first two phenotypes are associated with the interactions; C: Power when the first three phenotypes are associated with the interactions; D: Power when all four phenotypes are associated with the interactions.

Fig 5

Statistical power of Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI under moderate among-phenotype correlation ρ = 0.5.

The horizontal and vertical axes represent the proportion of causal rare variants and the statistical power, respectively. A: Power when only the first phenotype is associated with the interactions; B: Power when the first two phenotypes are associated with the interactions; C: Power when the first three phenotypes are associated with the interactions; D: Power when all four phenotypes are associated with the interactions.

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 exploit GEI pleiotropy effectively. For instance, in Fig 6, with the causal proportion θ = 0.2, when only the first phenotype is associated with the interactions, the powers of Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI are 0.013, 0.206, 0.292 and 0.045, respectively, as shown in Fig 6A. When the first two phenotypes are associated with the interactions, the corresponding power improves to 0.092, 0.607, 0.629 and 0.152, as shown in Fig 6B. With three phenotypes associated with the interactions, the power become even larger and are 0.173, 0.766, 0.767 and 0.252, as shown in Fig 6C. When all phenotypes are associated with the interactions, the power further increase to 0.236, 0.827, 0.828 and 0.345, as shown in Fig 6D.

We can also see from Figs 4–6 that Hom-GEI, Het-GEI and PPK-GEI have enhanced power as the correlation among phenotypes becomes stronger, but LPK-GEI suffers power loss. For instance, we observe the power of the four tests with causal proportion θ = 0.2. When the correlation among phenotypes is weak, the power values of Hom-GEI, Het-GEI, PPK-GEI and LPK-GEI are 0.105, 0.508, 0.523 and 0.370, respectively, in Fig 4C. When the correlation is moderate, the power values are 0.104, 0.597, 0.598 and 0.311, as shown in Fig 5C. With a strong correlation, the power values are 0.173, 0.766, 0.767 and 0.252, as shown in 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 single-phenotype 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 TOWGE-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.

Fig 7

Statistical power of Hom-GEI, Het-GEI, PPK-GEI, LPK-GEI and TOWGE-FCT with different correlated phenotypes, ρ = 0.5, θ = 0.1.

The horizontal and vertical axes represent the number of phenotypes associated with interactions and the statistical power, respectively. The significance level is 0.05.

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, age2, 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 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.

We carried out genome-wide analysis on 18,101 genes from 22 autosomal chromosomes. Variants in the genotype dataset were annotated via VEP [30]. We restricted to variants annotated as stop_loss, missense_variant, start_lose, splice_donor_variant, inframe_deletion, frameshift_variant, splice_acceptor_variant, stop_gained or inframe_insertion with PolyPhen scores larger than 0.15 and Sift scores less than 0.05 [31]. Those variants with MAFs less than 3% were extracted via PLINK [32]. Genotypes were further transformed into numeric values using fcGENE [33].

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 consider a suggestive significance level at 1×10−4, twelve genes passed the threshold, whose details are presented in Table 1.

Fig 8

Manhattan plot of genome-wide multiphenotype analysis of gene-Hb interactions in BPs.

The horizontal and vertical axes represent the genomic position and the negative log10 of the p-values, respectively. A: Het-GEI; B: PPK-GEI.

Fig 9

Q–Q plots of genome-wide multiphenotype GEI analysis of gene-Hb interactions in BPs.

The horizontal and vertical axes represent the negative log10 of the expected p-values and the negative log10 of the observed p-values, respectively. A: Het-GEI; B: PPK-GEI.

Table 1

Genes showing suggestive evidence (p-values < 1 × 10−4) of gene-Hb interactions in the Het-GEI or PPK-GEI tests.

Gene	Chr	Pos	RV Num	P value
				Het-GEI	PPK-GEI
ALX3	1	110060374: 110070700	27	3.34×10−5	5.21×10−5
CFHR1	1	196819730: 196832189	19	4.96×10−5	2.50×10−4
ARID5A	2	96536726: 96552634	16	9.60×10−2	1.53×10−4
OXCT1	5	41730064: 41870689	18	2.29×10−5	5.99×10−5
EMC2	8	108443623: 108486907	11	6.69×10−5	8.32×10−5
COQ4	9	128322511: 128334072	36	2.12×10−1	6.26×10−4
FAM120AOS	9	93446499: 93453592	5	3.57×10−5	3.01×10−6
NCR3LG1	11	17351761: 17377321	6	5.59×10−4	2.30×10−5
OR2D3	11	6921001: 6921994	15	1.36×10−1	3.86×10−4
CPSF6	12	69239536: 69274358	4	3.81×10−5	4.76×10−6
PYDC1	16	31215961: 31217074	2	5.29×10−5	2.41×10−4
LEUTX	19	39776593: 39786135	15	2.20×10−6	7.43×10−6

Abbreviations are as follows: Chr, chromosome. Pos, position. RV Num, number of rare variants.

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 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.

Table 2

Multiphenotypic analyses and single-phenotype analyses of gene-Hb interactions in BP-associated genes.

Gene	Chr	Pos	RV Num	P value
				Het-GEI	PPK-GEI	INT-FIX (SBP)	INT-FIX (DBP)
NPR1	1	153678687: 153693992	75	0.459	0.528	0.688	0.291
AGT	1	230702522: 230714590	46	0.556	0.497	0.599	0.213
PHC3	3	170087579: 170181749	40	0.630	0.371	0.264	0.157
NR3C2	4	148078763: 148442520	17	0.922	0.690	0.542	0.240
SLC39A8	4	102251040: 102345498	23	0.766	0.522	0.296	0.289
ENPEP	4	110476072: 110563337	74	0.663	0.785	0.510	0.912
NOS3	7	150991055: 151014599	126	0.183	0.065	0.129	0.026
ZFAT	8	134477787: 134713049	67	0.091	0.161	0.353	0.309
DBH	9	133636362: 133659344	81	0.828	0.627	0.371	0.339
PLCE1	10	93993988: 94328391	73	0.269	0.343	0.495	0.269
PLCB3	11	64251522: 64269452	85	0.227	0.344	0.241	0.792
ALKBH8	11	107502726: 107565735	41	0.924	0.926	0.644	0.871
PDE3B	11	14643722: 14872058	81	0.336	0.542	0.780	0.500
PDE3A	12	20369244: 20684107	26	0.205	0.304	0.440	0.400
RAPGEF3	12	47734669: 47759106	63	0.051	0.103	0.106	0.888
HDAC7	12	47782723: 47819980	38	0.402	0.470	0.673	0.336
HEATR5A	14	31291787: 31420582	129	0.740	0.468	0.165	0.382
SOS2	14	50117127: 50231381	60	0.170	0.491	0.647	0.981
SLC9A3R2	16	2026867: 2039026	47	0.968	0.893	0.685	0.594
MYO1C	17	1464185: 1492707	165	0.468	0.038	0.018	0.011
GATA5	20	62463496: 62475970	77	0.610	0.704	0.756	0.475
FOXS1	20	31844299: 31845617	44	0.563	0.553	0.285	0.592

Abbreviations are as follows: Chr, chromosome. RV Num, number of rare variants. Pos, position. SBP, systolic blood pressure. DBP, diastolic blood pressure.

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-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.

15 Jun 2022

6 Jul 2022

Respond to Editors:

We are herein resubmitting our revised manuscript titled “Kernel-Based Gene–Environment Interaction Tests for Rare Variants with Multiple Quantitative Phenotypes” (Manuscript Number: PONE-D-22-11234) to PLoS One.

We would like to thank you for giving us the opportunity to revise our manuscript and respond to the comments raised by the reviewer. We also want to thank the reviewer for his detailed comments and constructive suggestions, which help improving the quality of this paper substantially. The manuscript has been revised according to the comments of the reviewer and a point-by-point response to the reviewer’s comments is enclosed.

The resubmitted manuscript has been formatted according to the style of PLoS One. We have checked our reference list and added four references in the revised manuscript, the four references were cited in new contents that have added in the revised manuscript according to the reviewer’s comments. ORCID ID has been linked to our Editorial Manager account. Besides, we have added “Author Contributions” in our revised manuscript. Any changes made to our original manuscript have been highlighted in our revised manuscript. The “Funding” statement and “Data Availability Statement” are corrected as follows:

Funding

This work was supported by the national Thousand Youth Talents Plan to GS. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability Statement

This research has been conducted using the UK Biobank Resource under Application Number 44080. The genetic and phenotype datasets are available via the UK Biobank data access process. More details are available at http://www.ukbiobank.ac.uk/register-apply/.

Respond to Reviewer:

1.Have you considered comparing your four proposed methods with current existing methods in GEI pleiotropy field via simulations? Also, can you discuss what are the advantages of the proposed models and limitations of the existing methods?

Response: Thank you for the suggestion! We conducted an extensive literature search and found only one method available for testing GEI with rare variants and multiple phenotypes (Zhang et al. 2019). 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 gene-environment interactions for rare variants (TOW-GE) (Zhao et al. 2020); 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 TOWGE-FCT would become larger with the increasing number of phenotypes, which might limit statistical power of the test.

In the revised manuscript, we have added simulations and compared our proposed four multiphenotype GEI tests with TOWGE-FCT. The results showed that all the five tests have enhanced power as the number of phenotypes associated with interactions increases, our proposed Het-GEI, PPK-GEI and LPK-GEI had larger power than TOWGE-FCT, but Hom-GEI had smaller power than TOWGE-FCT. Moreover, TOWGE-FCT uses a permutation test to estimate p value for each phenotype, which is very time consuming. In comparison, our proposed four tests have asymptotic distributions and p values can be computed analytically. Thus, our proposed tests have much smaller computational cost. We also added discussion on the advantages of our proposed tests and limitations of TOWGE-FCT in the revised manuscript.

Zhang J, Sha Q, Hao H, Zhang S, Gao XR, Wang X. Test Gene-Environment Interactions for Multiple Traits in Sequencing Association Studies. Hum Hered. 2019;84(4-5):170-196. https://doi.org/10.1159/000506008 PMID: 32417835

Zhao Z, Zhang J, Sha Q, Hao H. Testing gene-environment interactions for rare and/or common variants in sequencing association studies. PLoS One. 2020;15(3):e0229217. https://doi.org/10.1371/journal.pone.0229217 PMID: 32155162

2.In the UK Biobank data application, there are several genes identified by your proposed models to have interactions. Are these genes confirmed in past literature to show correlations with BP phenotypes? Have you tried anyway to validate the results as well as your models?

Response: Thank you for the suggestion! In our work, we identified three genes that have interactions with Hb in BP phenotypes. For NOS3 and MYO1C genes, they have been reported to be correlated with BP phenotypes in (Surendran et al. 2020).

For LEUTX gene, in an analysis of 454,787 UK Biobank participants, LEUTX was not associated with the BP phenotypes at the nominal significance level of 0.05 (Backman et al. 2021). While our proposed tests identified LEUTX that shows interaction with Hb in BP phenotypes, the result lacks independent validation. Thus, it should be considered as being preliminary and further experiments are necessary. We acknowledged the limitation of the result in the revised manuscript.

Surendran P, Feofanova EV, Lahrouchi N, Ntalla I, Karthikeyan S, Cook J, et al. Discovery of rare variants associated with blood pressure regulation through meta-analysis of 1.3 million individuals. Nat Genet. 2020;52(12):1314-1332. https://doi.org/10.1038/s41588-020-00713-x PMID: 33230300

Backman JD, Li AH, Marcketta A, Sun D, Mbatchou J, Kessler MD, et al. Exome sequencing and analysis of 454,787 UK Biobank participants. Nature. 2021;599(7886):628-634. https://doi.org/10.1038/s41586-021-04103-z PMID: 34662886

Submitted filename: Response to Reviewers.doc

Click here for additional data file.

26 Sep 2022

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

PONE-D-22-11234R1

Dear Dr. Shi,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Mehdi Rahimi, Ph.D.

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: (No Response)

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

**********

30 Sep 2022

PONE-D-22-11234R1

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

Dear Dr. Shi:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Mehdi Rahimi

Academic Editor

PLOS ONE