Testing gene-environment interactions for rare and/or common variants in sequencing association studies

The risk of many complex diseases is determined by a complex interplay of genetic and environmental factors. Advanced next generation sequencing technology makes identification of gene-environment (GE) interactions for both common and rare variants possible. However, most existing methods focus on testing the main effects of common and/or rare genetic variants. There are limited methods developed to test the effects of GE interactions for rare variants only or rare and common variants simultaneously. In this study, we develop novel approaches to test the effects of GE interactions of rare and/or common risk, and/or protective variants in sequencing association studies. We propose two approaches: 1) testing the effects of an optimally weighted combination of GE interactions for rare variants (TOW-GE); 2) testing the effects of a weighted combination of GE interactions for both rare and common variants (variable weight TOW-GE, VW-TOW-GE). Extensive simulation studies based on the Genetic Analysis Workshop 17 data show that the type I error rates of the proposed methods are well controlled. Compared to the existing interaction sequence kernel association test (ISKAT), TOW-GE is more powerful when there are GE interactions’ effects for rare risk and/or protective variants; VW-TOW-GE is more powerful when there are GE interactions’ effects for both rare and common risk and protective variants. Both TOW-GE and VW-TOW-GE are robust to the directions of effects of causal GE interactions. We demonstrate the applications of TOW-GE and VW-TOW-GE using an imputed data from the COPDGene Study.

1 The etiology of many diseases is characterized by the interplay between genetic and 2 environment factors. For example, anthracyclines are one of the most effective classes of 3 chemotherapeutic agents currently available for cancer treatment. The therapeutic 4 potential of anthracyclines, however, is limited because of their strong dose-dependent 5 relation with progressive and irreversible cardiomyopathy leading to congestive heart 6 failure. Both gene hyaluronan synthase 3 ( HAS 3) and gene CUGBP Elav-like family 7 member 4 (CELF 4) modify the risk of anthracycline on the development of 8 anthracycline-related cardiomyopathy [1,2]. A genome-wide gene environment 9 interaction analysis indicates that gene EBF 1 plays together with stress associated 10 October 4, 2019 1/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint with cardiovascular disease. Additionally, gene EBF 1 not only shows gene-by-stress 11 interaction effect for hip circumference but also indicates gene-by-stress interaction 12 effects for waist circumference, body mass index (BMI), fasting glucose, type II diabetes, 13 and common carotid intimal medial thickness (CCIMT) [3]. 14 To date, most of the successful findings in gene-environment (GE) interactions are 15 for common genetic variants. There has been very limited success in findings for rare 16 variants’ GE interactions. This is often attributed to study design issues, such as sample 17 size or population heterogeneity [4]. Lack of statistical methodology on rare variants’ 18 GE also contributes to the limitations. 19 Rare variants, which are usually defined as genetic variants with minor allele 20 frequency (MAF) less than 5% (or 1%), may play an important role in studying the 21 etiology of complex human diseases. Numerous statistical methods have been developed 22 for testing the main effects of rare variants, such as the sequence kernel association test 23 (SKAT) [5], the combined multivariate and collapsing (CMC) method [6], the weighted 24 sum statistic (WSS) [7], and Testing the effect of an Optimally Weighted combination of 25 variants (TOW) [8]. 26 To our knowledge, limited methods have been developed for testing GE interactions 27 in sequencing association studies. Existing methods for assessing common variants by 28 environment interactions, such as the gene-environment interactions association test 29 (GESAT) [9] are less powerful when naively applied to rare variants [10]. To test rare 30 variants by environment interactions, [10] developed the interaction sequence kernel 31 association test (ISKAT) to assess the effects of rare variants by environment 32 interactions. As ISKAT considers the special weights Beta(MAF; 1, 25), the beta 33 distribution density function with parameters 1 and 25 evaluated at the sample MAF, 34 which is the recommended weight for ISKAT when there is no prior information, ISKAT 35 may lose power when the MAFs of causal variants are not in the range (0.01,0.035) [11]. 36 In this article, to test for rare and/or common variants and environment interactions 37 in sequencing association studies, we develop two novel methods: 1) Testing the 38 Optimally weighted combination of GE interactions for rare variants (TOW-GE); 2) 39 testing effects of weighted combination of GE interactions for both rare and common 40 variants (variable weight TOW-GE, refer to this statistic as VW-TOW-GE). Both 41 TOW-GE and VW-TOW-GE are robust to directions of effects of causal GE 42 interactions. We evaluate the performance of the proposed methods via simulation 43 studies and real data analysis using the imputed sequencing data from the COPDGene 44 Study. 45

46 Consider n unrelated individuals sequenced in a testing region with m genetic variants. 47 In the testing region, we are interested in testing the effects of p rare variants (p<m) by 48 environment interactions on a trait, which can be a quantitative or qualitative trait. For 49 ease of presentation, we only consider a single environmental factor. The method can be 50 easily extended to the case when there are multiple environmental factors. For 51 individuals i = 1,...,n , let yi denote the trait, Xi = (xi1,...,x iq)T denote the q 52 covariates,Gi = (gi1,...,g ip)T denote genotypes for the p rare variants in a genomic 53 region (a gene or a pathway) and Ei as the environmental factor. Let 54 Si = (Eigi1,...,E igip)T be a vector of variants by environment interaction terms for 55 the ith individual. 56 We use the generalized linear model (GLM) to model the relationship between the trait values yi and covariatesXi, genotypes Gi, environmental factor Ei and GE October 4, 2019 2/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint interactions,Si: g(E(yi|Xi,Gi,Ei)) =XT i α1 +Eiα2 +GT iα3 +ST i β (1) =˜XT i α +ST i β where g(·) is a canonical link function. Two commonly used models under the 57 generalized linear model framework are the linear model with the identity link for a 58 continuous or quantitative trait, and the logistic regression model with the Logit link for 59 a binary trait. α1,α 2,α3, and β are defined as q× 1 coefficient vector of covariate, the 60 coefficient of environmental factor, p× 1 coefficient vector of genotype and p× 1 61 coefficient vector of GE interactions for the ith individual and the trait, respectively. 62 Let ˜Xi = (XT i ,Ei,GT i )T andα = (α1,α 2,α3)T . Testing the association between the 63 trait and the rare variants by environment interactions is equivalent to testing the null 64 hypothesis H0 :β = 0. 65 We develop a score test by treating α as nuisance parameters and then adjust both 66 the trait value yi andSi for the covariatesXi, the genotypic score Gi, and the 67 environmental variableEi by applying linear regression. Denote ˜yi as the residual of yi 68 and ˜Si = (˜si1,..., ˜sip) as the residual of Si, regressed on ˜Xi. Then, the relationship 69 between ˜yi and ˜Si can be modeled by the GLM: 70 g(E(˜yi|˜Si)) =β∗ 0 + ˜ST i β∗ (2) TestingH0 :β = 0 in (1) is equivalent to testing H0 :β∗ = 0 in (2) (Sha et al., [8]). Here, we utilize a weight selection scheme proposed by Sha et al. [8] on our new model to test the effect of a weighted combination of GE, ˜si =∑p j=1wj˜sij. Following Sha et al. [12], we propose the following score test statistic under the generalized linear model: S(w1,...,w p) =n (∑n i=1(˜yi− ˜y)(˜si− ˜s) )2 ∑n i=1(˜yi− ˜y)2∑n i=1(˜si− ˜s)2 =n (∑p j=1wj ∑n i=1(˜yi− ˜y)(˜sij− ˜sj) )2 ∑n i=1(˜yi− ˜y)2∑n i=1(˜si− ˜s)2 Because GE interactions for rare variants are essentially independent, we have: n∑ i=1 (˜si− ˜s)2 = p∑ j=1 p∑ l=1 wjwl n∑ i=1 (˜sij− ˜sj)(˜sil− ˜sl) ≈ p∑ j=1 w2 j n∑ i=1 (˜sij− ˜sj)2 Thus, as a function of (w1,...,w p), the score test statistic S(w1,...,w p) reaches its 71 maximumS0(w0 1,...,w 0 p) =n∑n i=1(˜yi− ˜y)(˜s0 i− ˜s 0 )/∑n i=1(˜yi− ˜y)2 when 72 w0 j =∑n i=1(˜yi− ˜y)(˜sij− ˜sj)/∑n i=1(˜sij− ˜sj)2 and ˜s0 i =∑p j=1w0 j ˜sij. 73 Similarly, we define the statistic to Test the effect of the Optimally Weighted 74 combination of GE interactions (TOW-GE), ∑p j=1w0 j ˜sij, as: 75 TTOW −GE = n∑ i=1 (˜yi− ˜y)(˜s0 i− ˜s 0 ) (3) which is equivalent to S0(w0 1,...,w 0 p), where ∑n i=1(˜yi− ˜y)2 can be viewed as a constant 76 when we use a permutation test to evaluate p-values. 77 October 4, 2019 3/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint The optimal weight w0 j is equivalent to w0∗ j =ρ(˜y, ˜sj)/ √∑n i=1(˜sij− ˜sj)2, where 78 ρ(˜y, ˜sj) is the correlation coefficient between ˜y = (˜y1,..., ˜yn) and ˜sj = (˜s1j,..., ˜snj). 79 From the expression of w0∗ j , we can see that it is proportional to ρ(˜y, ˜sj) and thus w0 j 80 will put large weights to the GE interactions that have strong associations with the trait 81 and also adjust for the direction of the association. Simultaneously, w0∗ j is proportional 82 to 1/ √∑n i=1(˜sij− ˜sj)2 and w0 j will put large weights to GE interactions with small 83 variations which are common in GE interactions for rare variants. 84 TOW-GE focuses primarily on rare variants by environment interactions and it may 85 lose power because of the small weights on common variants by environment 86 interactions. Thus, to test the GE interactions’ effects of both rare and common 87 variants, we propose the following variable weight TOW-GE denoted as VW-TOW-GE. 88 We first divide GE interactions into two parts based on rare or common variants and 89 then we apply TOW-GE to the two parts separately. Let 90 Tλ =λ Tr√ var(Tr) + (1−λ) Tc√ var(Tc) where Tr and Tc denote the test statistics of 91 TOW-GE for GE interactions’ effects of rare and common variants, respectively. λ is a 92 tuning parameter. Denote pλ as the p-value of Tλ, and then the test statistic of 93 VW-TOW-GE is defined as TVW −TOW −GE = min0≤λ≤1pλ. In this study, we use a 94 simple grid search method to choose the tuning parameter λ and minimize the p-value. 95 Divide the interval [0, 1] into K subintervals of equal-length. Let λk =k/K for 96 k = 0, 1,...,K . Then, min 0≤λ≤1pλ = min0≤k≤Kpλk. 97 The p-value of TVW −TOW −GE can be evaluated by permutation tests following similar permutation tests for variable weight TOW (VW-TOW) proposed by [8]. Suppose that we perform B times of permutations. In each permutation, we randomly shuffle the trait values. Let T (b) r and T (b) c denote the values of Tr and Tc, respectively, based on the bth permuted data, where b = 0 represents the original data. Based on T (b) r and T (b) c (b = 0, 1, 2,...,B ), we can calculate T (b) λk for b = 0, 1, 2,...,B and k = 0, 1, 2,...,K , where var(Tr) and var(Tc) are estimated using T (b) r and T (b) c (b = 0, 1, 2,...,B ). Then, we transfer T (b) λk to p(b) λk by p(b) λk = # { T (d) λk :T (d) λk >T (b) λk ford = 0, 1, 2,...,B } B Let p(b) = min0≤k≤Kp(b) λk . Then, the p-value of TVW −TOW −GE is given by # { p(b): p(b)<p(0)forb = 0, 1, 2,...,B } B Simulation 98 We compared the performance of our proposed methods with the interaction sequence 99 kernel association test (ISKAT) [10], the modified WSS for testing the effects of GE 100 interactions [7] and the modified CMC method for testing the effects of GE 101 interactions [6]. In this study, the rank sum test used by WSS and the T 2 test used by 102 CMC were replaced with the score test based on residuals ˜yi and ˜sij. The empirical 103 Mini-Exome genotype data provided by the GAW17 is used for simulation studies. The 104 dataset contains genotypes of 697 unrelated individuals on 3,205 genes. Because gene 105 ELAVL 4 in GAW17 was used to simulate GE interaction’s effect on quantitative trait 106 October 4, 2019 4/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint Q1 which follows a normal distribution, we chose gene ELAVL 4 in our simulation 107 study. Gene ELAVL 4 has 10 variants, containing 8 rare variants and 2 common 108 variants. Rare variants in the simulation are defined with MAF < 0.05. 109 To evaluate type I error, we generate trait values independent of GE interactions (e.g. β1 = 0 and βc = 0) by using the model: Y1 = 0.5X1 + 0.5X2 +Eα1 +GTα2 +STβ1 +Scβc 1 +ϵ1 where ϵ1 follows a normal distribution with mean as 0 and variances as σ2 1 = 1; 110 α1 = 0.015;S is GE interactions for rare variants and Sc is GE interaction for a 111 common variant. We consider two covariates: a standard normal covariate X1 and a 112 binary covariateX2 with P (X2 = 1) = 0.5. The environmental factor E is assumed to 113 be continuous following a standard normal distribution. 114 For type I error evaluation, we consider two different cases: 1) testing the effects of 115 GE interactions for rare variants ; 2) testing the effects of GE interactions for both rare 116 and common variants. For each case, we consider two scenarios: (a) with main effect; 117 (b) without main effect in the model. When the main effects exist, we set the 118 magnitudes of vector α2 as 0.3 and the sign of each coefficient is random sampled from 119 (−1, 1). When main effects do not exist, we set α2 = 0. 120 For power comparisons, the phenotype is generated using similar settings to type I 121 evaluation except for existing GE interactions’ effects. We compare the power of 122 TOW-GE, ISKAT, WSS and CMC to test rare variant GE interactions’ effects 123 considering two scenarios: (a) including main effects, α2̸= 0 for rare variants; (b) no 124 main effects, α2 = 0 for rare variants. We vary the number of non-zero in the vector βi, 125 the proportion of non-zero in βi that are positive, and the magnitudes of the non-zero 126 βij. We set the magnitudes of the non-zero βij’s as|βij| =c, and increase c from 0.1 to 127 0.5. In each simulation scenario, p-values are estimated by 10,000 permutations and 128 1,000 replicated samples. 129 Simulation results 130 The empirical type I error rates are shown in Table 1 and Table 2. For 10,000 replicated 131 samples, the 95% confidence intervals for type I error rates of nominal levels as 0.05, 132 0.01 and 0.001 are (0.046, 0.054), (0.008, 0.012) and (0.0004, 0.0016), respectively. 133 When there are (a) main effects, e.g. α2̸= 0, TOW-GE, VW-TOW-GE, ISKAT and 134 WSS control type I error rates well and the burden test CMC tends to have very 135 conservative type I error rates (top panel of Table 1 and Table 2). When there are (b) 136 no main effects. e.g. α2 = 0, all methods can control type I error rates well (bottom 137 panel of Table 1 and Table 2). 138 The results for testing the effects of GE interactions of rare variants when including 139 main effect and no main effect are given in Figure 1 and Figure 2, respectively. In both 140 of these two scenarios, we consider the sample size as 2000 without a GE interaction of a 141 common variant. We do not apply VW-TOW-GE here because it is designed for existing 142 GE interactions’ effects of both common and rare variants. The top, middle, and 143 bottom panels in Figure 1 and 2 provide results for three cases, e.g. when there are 2, 6 144 and 8 non-zero βij’s, respectively. The left and right panels of Figure 1 and 2 present for 145 two cases, e.g. 50% of the βij are positive and 100% of the βij are positive, respectively. 146 For each plot, we vary c, the magnitudes of the non-zero βij. As shown in the four plots 147 for the case when 50% of the βij are positive, TOW-GE is more powerful than the other 148 three tests. For the case when 100% of the βij are positive, WSS is relatively more 149 powerful than TOW-GE since all the GxEs have the same direction of effects. 150 October 4, 2019 5/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint TOW-GE is more powerful than the other two tests. However, WSS is very sensitive to 151 the directions of effects due to aggregation of GE interactions directly. Among the four 152 tests (TOW-GE, ISKAT, WSS and CMC) in the two different cases, CMC is the least 153 powerful test. CMC loses power as it gives GE interactions of common variants large 154 weights, and thus GE interactions of common neutral variants will introduce large noise. 155 Power comparisons of the five tests (TOW-GE, VW-TOW-GE, ISKAT, WSS and 156 CMC) for testing GE interaction effects for both rare and common variants are given in 157 Figure 3. For each scenario in Figure 3, we vary c from 0.02 to 0.1 and set 50% of the 158 βij as positive. Simultaneously, we set the coefficient of a common variant by 159 environment interactionβc i as positive and the magnitudes of βc i as twice of βij which is 160 the coefficient of a rare variant by environment interaction. From Figure 3, we can see 161 that VW-TOW-GE is the most powerful test. CMC is the second most powerful test as 162 CMC puts large weights on GE interactions of common variants and gains power 163 increment when the GE interaction of a common variant plays an important role as the 164 causal effect. WSS is the least powerful test, which loses power because it puts very 165 small weight on the GE interaction of the common variant. 166 TOW-GE, VW-TOW-GE, and ISKAT can all be considered as quadratic statistics 167 which have reasonable power across a wide range of alternative hypothesis. The three 168

are robust to the different directions of the GE interaction effects. We perform 169 a further assement for the three methods. Figure 4 shows the results. When there are 170 causal effects of GE interactions for both common and rare variants, VW-TOW-GE 171 outperforms TOW-GE and ISKAT. TOW-GE is more powerful than ISKAT except 172 when the magnitude of the GE interactions is less than 0.04. 173 Real data analysis 174 Chronic obstructive pulmonary disease (COPD) is one of the most common lung 175 diseases characterized by long term poor airflow and is a major public health 176 problem [13]. It is a complex disease which is influenced by genetic factors, 177 environmental influences, and genotype-environment interactions. We have known that 178 cigarette smoking is the major environmental determinant of COPD [14]. Several genes 179 have been suggested to play a role in the presence of a gene-by-smoking interaction term. 180 Specifically, [15] reported that the 30-repeat allele of HMOX 1 was associated with 181 COPD in presence of a gene-by-smoking (pack-years) interaction term. [14] presented 182 that the GSTM 1 gene was associated with severe chronic bronchitis in heavy smokers. 183 The COPDGene Study is a multi-center genetic and epidemiologic investigation to 184 study COPD [16]. This study is sufficiently large and appropriately designed for 185 analysis of COPD. In this study, we consider more than 5,000 non-Hispanic Whites 186 (NHW) participants where the participants have completed a detailed protocol, 187 including questionnaires, pre- and post-bronchodilator spirometry, high-resolution CT 188 scanning of the chest, exercise capacity (assessed by six-minute walk distance), and 189 blood samples for genotyping. The participants were genotyped using the Illumina 190 OmniExpress platform. The genotype data have gone through standard quality-control 191 procedures for genome-wide association analysis detailed at http: 192 //www.copdgene.org/sites/default/files/GWAS_QC_Methodology_20121115.pdf. 193 We imputed the COPD genotype data using the EUR haplotypes from the 1000 194 Genome Project as references. 195 Based on the literature of COPD [17,18], we selected 7 key quantitative 196 COPD-related phenotypes, including FEV1 (% predicted FEV1), Emphysema (Emph), 197 Emphysema Distribution (EmphDist), Gas Trapping (GasTrap), Airway Wall Area 198 (Pi10), Exacerbation frequency (ExacerFreq), Six-minute walk distance (6MWD), and 199 one qualitative phenotypes (case-control disease status denoted as COPD in following 200 October 4, 2019 6/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint tables). 3 covariates, including BMI, Age and Sex and one environmental factor 201 (Pack-Years) were considered in our analysis. EmphDist is the ratio of emphysema at 202 -950 HU in the upper 1/3 of lung fields compared to the lower 1 /3 of lung fields where 203 we did a log transformation on EmphDist in the following analysis, referred to [17]. In 204 the analysis, participants with missing data in any of these phenotypes were excluded. 205 To evaluate the performance of our proposed method on a real data set, we applied 206 all of the 5 methods (TOW-GE, ISKAT, WSS ,CMC, and VW-TOW-GE) to two 207 COPD associated genes (HMOX 1 and GSTM 1) through an interaction with cigarette 208 smoking. In the analysis, we removed the extreme rare SNPs (MAF <0.001) in any 209 genotypic variants and missing value in any of the 7 phenotypes and 3 covariates. We 210 considered three different scenarios: (1) main effect; (2) gene-by-smoking interaction 211 with main effect and (3) gene-by-smoking interaction without main effect. When we 212 considered only the main effect, we used five existing methods (TOW-GE, SKAT, WSS 213 ,CMC, and VW-TOW) which are specifically designed for testing the main effect of a 214 gene. We adopted 10 4 permutations for our methods and used 0.05 as the significance 215 level. The results for testing association between COPD and gene HMOX 1 and 216 GSTM 1 are summarized in Table 3. At gene HMOX 1, both TOW-GE and modified 217 WSS verified significant GE intecation effects and with main effect for two traits Emph 218 and Pi10, ISKAT and VW-TOW-GE verified significant GE intecation effects and 219 without main effect for trait Emph. 220

221 Recent evidence shows that gene-environment interactions of rare variants may play an 222 important role in explaining the etiology of a complex disease. However, there are 223 limited methods that can be employed to test the effects of GE interactions for rare 224 variants. In this study, we propose two new methods for testing GE interactions for rare 225 variants only or for both rare and common variants. We employ a generalized linear 226 model to model the relationship between the trait and the GE interactions. Our model 227 focuses on GE interactions by first adjusting for genetic main effects, enviroumental 228 main effects, and possible covariates. Two methods are designed for different scenarios 229 through specific weigh-selection mechanisms. TOW-GE assigns the majority of weights 230 on rare variants by environment interactions. VW-TOW-GE balances common and rare 231 variants by performing weight assignments separately for common variants by 232 environment interactions and rare variants by environment interactions. Both methods 233 achieve the best possible power with an adaptive weight selection procedure. 234 In the application, we have tested genetic association for 7 traits of COPD. Our 235 proposed methods verified the most significant GE interactions, especially for 236 gene-by-smoking interactions without main effect and performs the best compared to 237 other methods. In simulation studies, we also demonstrated that our proposed methods 238 perform better in different scenario: with main effect and without main effect. Our 239

show that the proposed methods TOW-GE or VW-TOW-GE demonstrate better 240 power in most cases compared with competing methods. 241 The power of a test varies according to the number of GE interactions of rare or 242 common variants, the effect directions of GE interactions, and the MAFs of variants. 243 When substantial of GE interactions have opposite directions of effects, the quadratic 244 statistics TOW-GE, VW-TOW-GE, and ISKAT are powerful. When effects of GE 245 interactions of common variants play a primary role, CMC is more powerful than 246 ISKAT, WSS, and has similar power to VW-TOW-GE. 247 In our proposed method, the optimal weights of TOW-GE are derived analytically; 248 thus the computations cost is relatively small. On the other hand, TOW-GE is flexible 249 and allows for prior biological information to be incorporated by using flexible weights, 250 October 4, 2019 7/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint such as weights derived from the expression quantitative trait locus (eQTL), which may 251 further improve the power of TOW-GE. In addition, TOW-GE allows for adjustment of 252 covariates. The covariates could be demographic variables, environmental variables, 253 clinical variables, and/or principle components of genotype scores. The adjustment of 254 covariates makes TOW-GE not only able to eliminate the effect of confounders but also 255 able to correct for possible population stratification in admixed populations. One 256 possible advantage of TOW-GE compared to ISKAT is that TOW-GE utilizes the 257 residuals of both the trait value and the GE interactions, which are obtained by 258 adjusting for covariats from linear regression models, respectively, while ISKAT utilizes 259 only the residual of the trait value. 260

261 The Genetic Analysis workshops are supported by NIH grant R01 GM031575 from the 262 National Institute of General Medical Sciences. Preparation of the Genetic Analysis 263 Workshop 17 Simulated Exome Data Set was supported in part by NIH R01 MH059490 264 and used sequencing data from the 1000 Genomes Project (www.1000genomes.org). 265 This research used data generated by the COPDGene study (phs000179/HMB and 266 phs000179/DS-CS-RD), which was supported by National Institutes of Health (NIH) 267 grants U01HL089856 and U01HL089897. The content is solely the responsibility of the 268 authors and does not necessarily represent the official views of the National Heart, 269 Lung, and Blood Institute or the National Institutes of Health. The COPDGene project 270 is also supported by the COPD Foundation through contributions made by an Industry 271 Advisory Board comprised of Pfizer, AstraZeneca, Boehringer Ingelheim, Novartis, and 272 Sunovion. 273 A superior high-performance computing infrastructure at University of North Texas 274 was used in obtaining results presented in this publication. 275 Author Contributions 276 Formal analysis: Zihan Zhao, Jianjun Zhang 277 Methodology: Zihan Zhao, Han Hao 278 Visualization: Qiuying Sha 279 Writing-original draft: Zihan Zhao, Han Hao 280 Writing-review & editing: Zihan Zhao, Jianjun Zhang, Qiuying Sha, Han Hao 281 Figure legends 282 Fig 1. Power comparisons of the four tests (TOW-GE, ISKAT, WSS and CMC) for 283 testing GE interaction effects for rare variants on a continuous outcome when there are 284 main effects (n=2000 and the significance level of α = 0.05). 285 Fig 2. Power comparisons of the four tests (TOW-GE, ISKAT, WSS and CMC) for 286 testing GE interaction effects of rare variants on a continuous outcome when there are 287 no main effects (n=2000, significance level of α = 0.05). 288 Fig 3. Power comparisons of the five tests (TOW-GE, ISKAT, WSS,CMC and 289 VW-TOW-GE) for testing GE interaction effects for both rare and common variants on 290 a continuous outcome (n=2000 and the significance level of α = 0.05). Left panel: With 291 main effect; Right panel: With no main effect. 292 Fig 4. Power comparisons of the three quadratic tests (TOW-SE, iSKAT, and 293 VW-TOW-SE) for testing GE interaction effects of both rare and common variants on a 294 October 4, 2019 8/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint continuous outcome (n=2000, the significance of α = 0.05). Left panel: With main 295 effect; Right panel: Without main effect. 296 October 4, 2019 9/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint T able 1. Type 1 error rates for testing the effects of GE interactions of rare variants in the presence of main effects (top panel) and in the absence of main effects (bottom panel) (n=2000). With main effect α-level TOW-GE ISKAT WSS CMC n = 2000 0.05 0.042 0.066 0.047 0.027 0.01 0.01 0.017 0.011 0.005 0.001 0.001 0.009 0.000 0.000 Without main effect n = 2000 0.05 0.054 0.066 0.050 0.043 0.01 0.006 0.014 0.012 0.012 0.001 0.000 0.004 0.000) 0.000 T able 2. Type 1 error rates for testing the effects of GE interactions for both rare and common variants in the presence of main effects (top panel) and in the absence of main effects (bottom panel) (n=2000). With main effect α-level TOW-GE ISKAT WSS CMC VW-TOW-GE n = 2000 0.05 0.053 0.062 0.055 0.040 0.052 0.01 0.007 0.013 0.017 0.009 0.011 0.001 0.002 0.002 0.001 0.000 0.002 Without main effect n = 2000 0.05 0.051 0.056 0.048 0.049 0.058 0.01 0.006 0.012 0.012 0.011 0.014 0.001 0.001 0.003 0.001 0.001 0.002 October 4, 2019 10/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint T able 3. Summary results of association analysis for HMOX 1 based on the COPD dataset. The p-values are shown for testing the gene’s main effect (top panel), gene-by-smoking interaction with main effect (middle panel), gene-by-smoking interaction without main effect (bottom panel). Gene’s main effect trait TOW SKAT WSS CMC VW-TOW GasT rap 0.3917 0.542 0.2618 0.9038 0.5539 ExacerFreq 0.7050 0.6149 0.1922 0.9845 0.8155 Emph 0.0204 0.0328 0.0036 0.9155 0.0360 Pi10 0.0283 0.0266 0.0457 0.4501 0.0559 EmphDist 0.8083 0.7363 0.5172 0.7520 0.8888 6MWD 0.8299 0.8451 0.8526 0.9985 0.8642 FEV1 0.6825 0.6928 0.7057 0.6906 0.7691 COPD 0.8637 0.8277 0.8345 0.7540 0.8677 Gene-b y-smoking interaction with main effect trait TOW-GE ISKAT WSS CMC VW-TOW-GE GasT rap 0.7432 0.8001 0.2610 0.8894 0.8033 ExacerFreq 0.5883 0.2389 0.2768 0.9964 0.3921 Emph 0.4024 0.2718 0.1140 0.9861 0.5696 Pi10 0.1208 0.4084 0.0821 0.9948 0.0651 EmphDist 0.5315 0.4794 0.6006 0.9892 0.4886 6MWD 0.6174 0.3624 0.4211 0.9929 0.6793 FEV1 0.8656 0.7748 0.4419 0.9575 0.9178 COPD 0.2302 0.3029 0.9089 0.9424 0.3394 Gene-b y-smoking interaction without main effect trait TOW-GE ISKAT WSS CMC VW-TOW-GE GasT rap 0.3388 0.5724 0.1207 0.6040 0.4967 ExacerFreq 0.3818 0.2810 0.0915 0.9320 0.4513 Emph 0.0189 0.0487 0.0011 0.8288 0.0349 Pi10 0.0304 0.0532 0.0118 0.5587 0.0571 EmphDist 0.8166 0.8062 0.7610 0.8066 0.8217 6MWD 0.7253 0.3463 0.6810 0.9929 0.6811 FEV1 0.5604 0.7387 0.4519 0.3043 0.7280 COPD 0.8869 0.8877 0.8657 0.3204 0.9300 Note: The bold numbers represent p-values of significant tests (significance level = 0.05). October 4, 2019 11/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint T able 4. Summary results of association analysis for GSTM1 based on the COPD dataset. The p-values are shown for testing the gene’s main effect (top panel), gene-by-smoking interaction with main effect (middle panel), gene-by-smoking interaction without main effect (bottom panel) Gene’s main effect trait TOW SKAT WSS CMC VW-TOW GasT rap 0.2309 0.6152 0.6163 0.2479 0.2848 ExacerFreq 0.7198 0.7823 0.7677 0.9138 0.6594 Emph 0.1401 0.3901 0.8496 0.1686 0.2355 Pi10 0.0177 0.1069 0.1705 0.0749 0.0151 EmphDist 0.0077 0.0256 0.3545 0.0487 0.0082 6MWD 0.6011 0.6856 0.8401 0.9260 0.5707 FEV1 0.1190 0.5920 0.9013 0.2301 0.0866 COPD 0.2144 0.2178 0.4699 0.2078 0.3243 Gene-b y-smoking interaction with main effect trait TOW-GE ISKAT WSS CMC VW-TOW-GE GasT rap 0.8652 0.7482 0.5358 0.1096 0.9158 ExacerFreq 0.7417 0.0867 0.3860 0.0599 0.5606 Emph 0.6829 0.9901 0.6833 0.2927 0.7207 Pi10 0.2757 0.5465 0.4808 0.1164 0.4506 EmphDist 0.1287 0.2314 0.6639 0.6781 0.1126 6MWD 0.8144 0.8769 0.8893 0.3781 0.8384 FEV1 0.9389 0.4145 0.6640 0.1277 0.9169 COPD 0.9944 0.8870 0.7842 0.2098 0.9878 Gene-b y-smoking interaction without main effect trait TOW-GE ISKAT WSS CMC VW-TOW-GE GasT rap 0.5160 0.2723 0.8413 0.6725 0.5769 ExacerFreq 0.5887 0.7348 0.6194 0.2701 0.6796 Emph 0.1041 0.1114 0.6514 0.4787 0.1691 Pi10 0.0697 0.1112 0.1282 0.0844 0.0631 EmphDist 0.0071 0.0229 0.6162 0.1078 0.0131 6MWD 0.7759 0.9342 0.9903 0.8683 0.7867 FEV1 0.2833 0.4709 0.6934 0.2673 0.2254 COPD 0.3641 0.1615 0.4693 0.5593 0.4934 Note: The bold numbers represent p-values of significant tests (significance level = 0.05). October 4, 2019 12/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint

1. Wang X, Liu W, Sun C, Armenian S H, Hakonarson H, Hageman L, et al. Hyaluronan synthase 3 variant and anthracycline-related cardiomyopathy: a report from the children’s oncology group. Journal of Clinical Oncology. 2014;32(7):647–653. 2. Wang X, Sun C, Quiones-Lombraa A, Singh P, Landier W, Hageman L, et al. CELF4 variant and anthracycline-related cardiomyopathy: a children’s oncology group genome-wide association study. Journal of Clinical Oncology. 2016;34(8):863–870. 3. Singh A, Babyak MA, Nolan D K, Brummett B H, Jiang R, Siegler I C, et al. Gene by stress genome-wide interaction analysis and path analysis identify EBF1 as a cardiovascular and metabolic risk gene. European Journal of Human Genetics. 2015;23(6):854–862. 4. Thomas D. Gene-environment-wide association studies: emerging approaches. Nature Reviews Genetics. 2010;11(4):259-272. 5. Wu M C, Lee S, Cai T, Li Y, Boehnke M, Lin X. Rare-variant association testing for sequencing data with the sequence kernel association test. The American Journal of Human Genetics. 2011;89(1):82-93. 6. Li B, Leal S M. Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data. The American Journal of Human Genetics. 2008;83(3):311–321. 7. Madsen B E, Browning S R. A groupwise association test for rare mutations using a weighted sum statistic. PLoS genetics. 2009;5(2),e1000384. 8. Sha Q, Wang X, Wang X, Zhang S. Detecting association of rare and common variants by testing an optimally weighted combination of variants. Genetic epidemiology. 2012;36(6):561–571. 9. Lin X, Lee S, Christiani D C, Lin X. Test for interactions between a genetic marker set and environment in generalized linear models. Biostatistics. 2013;14(4):667–681. 10. Lin X, Lee S, Wu M C, Wang C, Chen H, Li Z, et al. Test for rare variants by environment interactions in sequencing association studies. Biometrics. 2016;72(1),156–164. 11. Yang X, Wang S, Zhang S, Sha Q. Detecting association of rare and common variants based on cross-validation prediction error. Genetic epidemiology. 2017;41(3):233–243. 12. Sha Q, Zhang Z, Zhang S. An improved score test for genetic association studies. Genetic epidemiology. 2011;35(5):350–359. 13. Murphy T F, Sethi S. Chronic obstructive pulmonary disease. Aging. 2002;19(10):761–775. 14. Sandford A J, Silverman E K. Chronic obstructive pulmonary disease 1: Susceptibility factors for COPD the genotype environment interaction. Thorax. 2002;57(8),763-741. October 4, 2019 13/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint 15. Hersh C P, DeMeo D L, Lange C, Litonjua A A, Reilly J J, Kwiatkowski D, et al. Attempted replication of reported chronic obstructive pulmonary disease candidate gene associations. American journal of respiratory cell and molecular biology. 2005;33(1),71–78. 16. Regan E A, Hokanson J E, Murphy J R, Make B, Lynch D A, Beaty T H, et al. Genetic epidemiology of COPD (COPDGene) study design. COPD: Journal of Chronic Obstructive Pulmonary Disease. 2011;7(1):32–43. 17. Chu J H, Hersh C P, Castaldi P J, Cho M H, Raby B A, Laird N, et al. Analyzing networks of phenotypes in complex diseases: methodology and applications in COPD. BMC systems biology. 2014;8(1):78. 18. Han M K, Kazerooni E A, Lynch D A, Liu L X, Murray S, Curtis J L, et al. Chronic obstructive pulmonary disease exacerbations in the COPDGene study: associated radiologic phenotypes. Radiology. 2011;26(1):274–282. October 4, 2019 14/14 .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint .CC-BY-NC-ND 4.0 International licensea certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under The copyright holder for this preprint (which was notthis version posted October 7, 2019. ; https://doi.org/10.1101/796540doi: bioRxiv preprint