Variance-component-based meta-analysis of gene–environment interactions for rare variants

Abstract Complex diseases are often caused by interplay between genetic and environmental factors. Existing gene–environment interaction (G × E) tests for rare variants largely focus on detecting gene-based G × E effects in a single study; thus, their statistical power is limited by the sample size of the study. Meta-analysis methods that synthesize summary statistics of G × E effects from multiple studies for rare variants are still limited. Based on variance component models, we propose four meta-analysis methods of testing G × E effects for rare variants: HOM-INT-FIX, HET-INT-FIX, HOM-INT-RAN, and HET-INT-RAN. Our methods consider homogeneous or heterogeneous G × E effects across studies and treat the main genetic effect as either fixed or random. Through simulations, we show that the empirical distributions of the four meta-statistics under the null hypothesis align with their expected theoretical distributions. When the interaction effect is homogeneous across studies, HOM-INT-FIX and HOM-INT-RAN have as much statistical power as a pooled analysis conducted on a single interaction test with individual-level data from all studies. When the interaction effect is heterogeneous across studies, HET-INT-FIX and HET-INT-RAN provide higher power than pooled analysis. Our methods are further validated via testing 12 candidate gene–age interactions in blood pressure traits using whole-exome sequencing data from UK Biobank.


Introduction
With the advancement of high-throughput sequencing technologies (Ansorge 2009), a variety of statistical tests (Li and Leal 2008;Madsen and Browning 2009;Wu et al. 2011;Lee et al. 2012;Cheng et al. 2014;Sun et al. 2016;Liu et al. 2019) have been developed to identify associations between rare variants and complex diseases or traits. For most complex diseases or traits, thousands of genetic variants identified by genome-wide association studies (GWASs) explain only a small proportion of their heritabilities, leaving much of those heritabilities still unexplained; this phenomenon is known as "missing heritability" (Tong and Hernandez 2020). Variants identified through GWASs are mostly common variants with minor allele frequencies (MAFs) larger than 5%. Studies have demonstrated that rare variants may explain some of the "missing heritability" Kao et al. 2017;Yu et al. 2018). For example, two rare variants have been identified to be associated with low-density lipoprotein cholesterol and high-density lipoprotein cholesterol (Igartua et al. 2017). Five rare variants have been associated with low systolic blood pressure (SBP) (He et al. 2017). On the other hand, complex diseases are often caused by a combination of genetic factors, environmental exposures, and the interplay between them; thus, gene-environment interaction (G Â E) may explain some of the "missing heritability" of complex diseases (Lim et al. 2020).
Recently, an increasing number of statistical methods have been developed for testing G Â E with rare variants. These methods aim to detect individually weak but collectively strong G Â E effects in a functional set, typically a gene or a gene region. Tzeng et al. (2011) proposed a similarity-based regression approach (SIMreg) to test the G Â E of rare variants in which trait similarity is regressed on pairwise genetic similarity. The approach is applicable for both continuous and binary traits and allows testing for main genetic effects, interaction effects, and joint effects. Jiao et al. (2013) presented a gene-based G Â E test (SBERIA) for casecontrol studies in which gene-environment correlations are used to filter out variants showing no promising G Â E effects. Lin et al. (2013) proposed a G Â E set association test (GESAT) under generalized linear models and tested the interaction effect between a marker set and an environmental variable for continuous and discrete traits. GESAT assumes the G Â E effect to be random and employs a variance component score test in a linear mixed model framework. Chen et al. (2014) proposed two G Â E tests for rare variants: INT-FIX, which treats main genetic effects as fixed effects, and INT-RAN, which treats main genetic effects as random. Lin et al. (2016) proposed the interaction sequence kernel association test (iSKAT) for assessing rare variants by environmental interactions, which is robust to the proportion of causal variants in a gene and the directions of the variants' interaction effects. Liu et al. (2016) developed a unified gene-based G Â E test, which filters out variants with little evidence of interaction effects. This method was shown to improve the power of interaction analysis by means of an optimal filtering threshold.
All these methods are designed for the analysis of interactions in a single study, which is often underpowered due to limited sample size. A very large number of samples are indeed necessary for identifying genetic variants, especially for finding G Â Es (Smith and Day 1984;McCarthy et al. 2008). To overcome this limitation, the sharing of results among consortia on the same disease or trait and meta-analyses that combine the results of multiple studies are routine practices (Panagiotou et al. 2013). A meta-analysis that combines summary results rather than individual-level raw data from multiple independent studies offers an increased effective sample size and boosted statistical power (Evangelou and Ioannidis 2013;Hu et al. 2013;Shi and Nehorai 2017;Jin and Shi 2019a, b). However, little work has been done on the meta-analysis of G Â Es with rare variants. Wang et al. (2018) suggested a gene-based meta-analysis with filtering to detect G Â E effects in case-control studies. In this approach, a variant-level "meta-filtering" test is first conducted, and metaanalysis techniques are then applied to test gene-based G Â E effects on the retained variants.
INT-FIX and INT-RAN are gene-based G Â E tests for rare variants in a single study, which were proposed by Chen et al. (2014) and are implemented in the rareGE package. In this study, we develop corresponding meta-analysis methods based on INT-FIX and INT-RAN, which are two powerful and accessible methods for testing G Â Es at the study level. Our meta-analysis methods combine gene-level summary statistics from studies via INT-FIX or INT-RAN tests and consider the G Â E effect to be either homogeneous or heterogeneous across studies. In simulation studies, we have evaluated our methods by examining their null distributions and statistical power. Our methods have been further applied in testing 12 candidate gene-age interactions in blood pressure (BP) regulation (Simino et al. 2014) using whole-exome sequencing data from UK Biobank (Zhao et al. 2020).

INT-FIX and INT-RAN in individual studies
Suppose that there are K studies in a meta-analysis for testing gene-based G Â E effects of rare variants in a region. For the k-th study, n k individuals have been sequenced in the region, which has m k variants, and y ki and G ki ¼ ðg ki1 ; g ki2 ; Á Á Á ; g kim k Þ denote the phenotype and genotypes of individual i ð1˚i˚n k Þ. Under the additive genetic model, g kij ¼ 0, 1, or 2 are the numbers of minor alleles. The first element of the covariate vector X ki ¼ ðX ki1 ; X ki2 ; Á Á Á ; X ki;ðp k þ1Þ Þ is the intercept, with a value of 1, and the other p k elements are covariates. The environmental factor E ki is included as one of the covariates for adjusting its main effect.
Under the assumption of independent observations, consider the following linear mixed model for testing the gene by E ki interaction as follows: where e ki $ Nð0; r 2 k Þ is an error term and a k ¼ ða k1 ; a k2 ; Á Á Á ; a kðp k þ1Þ Þ T represents the effects of the intercept and covariates. b k ¼ ðb k1 ; b k2 ; . . . ; b km k Þ T and c k ¼ ðc k1 ; c k2 ; . . . ; c km k Þ T are the main genetic effects and G Â E effects, respectively. Assume that c k has a mean of 0 and a covariance matrix sI m k , where I m k denotes the m k Â m k identity matrix. W k1 and W k2 denote the m k Â m k diagonal weight matrices with elements w kj1 and w kj2 (1˚j˚m k ), which are the weights of the main genetic effects and G Â E effects, respectively.
Let y k ¼ ðy k1 ; y k2 ; . . . ; y kn k Þ T be the phenotype vector, let E k ¼ diagfE ki g be an n k Â n k diagonal matrix, and let e k ¼ ðe k1 ;e k2 ; Á Á Á ;e kn k Þ T be the error vector. The linear mixed model is written in matrix form as The null hypothesis with no interaction effect for any variant in the gene is H 0 : s ¼ 0, and the alternative hypothesis of at least one variant having a nonzero interaction effect is H 1 : s > 0.
In the INT-FIX method (Chen et al. 2014), the main genetic effects are assumed to be fixed effects, and the b k are the regression coefficients of the main genetic effects. Under the null hypothesis, where Varðe k Þ ¼ V ¼ diagfr 2 k g. By means of maximum likelihood or restricted maximum likelihood functions, a k and b k can be estimated; the corresponding estimates are denoted by a^k F and b^k F , respectively. Then, the estimated mean and covariance matrix of y k are The INT-FIX statistics of testing the G Â E interaction is is the G Â E score statistic for the j-th variant. The asymptotic distribution of Q k;INTÀFIX is a mixture of chi-square distributions, and the P-value can be calculated via Kuonen's saddlepoint method (Kuonen 1999;Chen et al. 2014).
In the INT-RAN method (Chen et al. 2014), the main genetic effects are assumed to be random, and the elements of b k follow a normal distribution with zero mean and covariance matrix dI m k . Under the null hypothesis, the working model is written as (3), with Varðy k Þ ¼ dG k W k1 W k1 G T k þ r 2 k I n k , and a^k R , d^R, and r^k R 2 can be numerically estimated by means of maximum likelihood or restricted maximum likelihood functions (Chen et al. 2014). The estimated mean and covariance matrix of y k are is the G Â E score statistic for the j-th variant. The Q k;INTÀRAN statistic asymptotically follows a mixture of chi-square distributions, and the corresponding P-value can again be computed via Kuonen's saddlepoint method (Kuonen 1999;Chen et al. 2014).

Meta-analysis of summary statistics from INT-FIX and INT-RAN
The meta-analysis of G Â E effects is performed based on summary statistics from each study. The required summary statistics for each study include the MAFs of all variants, the G Â E score statistic (5) or (7) for each variant and the regional G Â E relationship matrix (Zhang and Lin 2003;Wu et al. 2011).
For the meta-analysis results of INT-FIX, the regional G Â E relationship matrix of the k-th study is For INT-RAN, the regional G Â E relationship matrix of the k-th study is For convenience, we assume that all variants can be observed in each of the K studies, that is, m 1 ¼m 2 ¼ Á Á Á ¼m K ¼m. When some variant is monomorphic in a study, its G Â E score statistic and corresponding elements in the relationship matrices can be simply set to zero.
We consider two types of meta-analyses, under scenarios in which the G Â E effects are homogeneous and heterogeneous across different studies. For the scenario in which the G Â E effects are homogeneous across the different studies, we propose the following G Â E test statistics for the meta-analysis of summary statistics from INT-FIX and INT-RAN: Similar to the meta-analysis of homogeneous main genetic effects presented in Lee et al. (2013), these two statistics first combine the weighted G Â E score statistics across different studies for each variant and then aggregate the squared score statistics of all variants in a gene. Here, the weights are based on the MAFs of the variants estimated across all studies. Under the null hypothesis, Q HOMÀINTÀFIX asymptotically follows a mixed 1-df chisquare distribution Chen et al. 2014), namely, Q HOMÀINTÀFIX $ P m j¼1 k j v 2 j;1 , where k 1 ;k 2 ; Á Á Á ;k m are the nonzero eigenvalues of and the P-value can be calculated via Kuonen's saddlepoint method (Kuonen 1999). Similarly, we have the test statistic Q HOMÀINTÀRAN $ P m j¼1 k j v 2 j;1 , where k 1 ;k 2 ; Á Á Á ;k m are the nonzero eigenvalues of For the scenario in which the G Â E effects are assumed to be heterogeneous across the different studies, we propose the following G Â E test statistics for the meta-analysis of the summary statistics from INT-FIX and INT-RAN: As in the meta-analysis of heterogeneous main genetic effects presented in Lee et al. (2013), these test statistics aggregate the squared and weighted score statistics across all studies and all variants. In HET-INT-FIX and HET-INT-RAN, the weights are based on the study-specific MAFs. Here, Q HETÀINTÀFIX $ P m j¼1 k j v 2 j;1 , where k 1 ;k 2 ; Á Á Á ;k m are the nonzero eigenvalues of W INTÀFIX , and Q HETÀINTÀRAN $ P m j¼1 k j v 2 j;1 , where k 1 ;k 2 ; Á Á Á ;k m are the nonzero eigenvalues of W INTÀRAN . The P-values of the mixed v 2 s can also be calculated via Kuonen's saddlepoint method (Kuonen 1999).
For trans-ethnic meta-analyses, we assume that the G Â E effects in studies of the same ancestry are homogeneous, and heterogeneous for studies of different ancestries. Suppose that there are K studies from B ancestries. In the b-th ancestry, there are as the number of studies from ancestry 1 to b, in particular, K~0 ¼ 0. We propose the following G Â E test statistics for the trans-ethnic meta-analysis of the summary statistics from INT-FIX and INT-RAN: The test statistics combine the weighted score statistic of the j-th variant from studies of the same ancestry, then aggregate the squared scores across all ancestries and all variants. Here, the weights are based on ancestry-specific MAFs. When all studies in the meta-analysis are of the same ancestry, Equations (16) and (17) reduce to Equations (10) and (11), respectively. When all studies have different ancestries, Equations (16) and (17) reduce to Equations (14) and (15), respectively.

Data availability
This research has been conducted using the UK Biobank Resource under Application Number 44080. Corresponding R codes for meta-analysis methods of testing G Â E effects in this article are available at GitHub: https://github.com/jlyx53/Codesfor-Meta-analyses-of-G-E-effects.

Numerical simulations
We conducted simulation studies to evaluate the null distributions and statistical power of our four meta-statistics for a continuous phenotype. We considered three scenarios as summarized in Table 1. Each has three studies with sample sizes of 1600, 2200, and 3200, referred to study 1, study 2, and study 3, respectively. For scenario 1, all studies are European (EUR) samples and have the same set of covariates; For scenario 2, three studies are all EUR samples but have different set of covariates; For scenario 3, three studies have different set of covariates and have different ancestries: the first two studies are EUR samples and the third study is African-American (AA) samples. In the power analysis, for each scenario, we considered cases with homogeneous and heterogeneous G Â E effects across studies, and also five different levels of main and interaction effects.
Using the calibrated coalescent model implemented in COSI (Schaffner et al. 2005), we first simulated 10,000 EUR haplotypes and 10,000 AA haplotypes over a 200 kb region. The simulation parameters were set to mimic the linkage disequilibrium structure, local recombination rate and population history of the EUR and AA populations ). Then, for scenarios 1 and 2, we randomly paired the EUR haplotypes to obtain diploid genotype data of 10,000 EUR individuals and randomly selected 7000 out of the 10,000 individuals. For scenario 3, we randomly paired the EUR haplotypes to obtain diploid genotype data of 10,000 individuals and randomly selected 3800 out of them. Similarly, we obtained diploid genotype data of 10,000 AA individuals, out of which we randomly selected 3200 individuals. Since the average exon length of a gene is approximately 3 kb (Pruitt et al. 2012), we randomly selected a 3 kb subregion within the 200 kb region for each replicate of the genotype data and retained variants with MAFs less than 0.03. We repeated this process 1000 times and generated 1000 replicates of genotype data sets for the three scenarios. On average, there were 54 rare variants in the 3 kb region under scenarios 1 and 2 and 87 under scenario 3. For scenarios 1 and 2, out of the genotype data of the 7000 EUR individuals, we randomly selected 1600, 2200, and 3200 for the samples of study 1, study 2, and study 3, respectively. For scenario 3, out of the genotype data of the 3800 EUR individuals, we randomly selected 1600 and 2200 for the samples of study 1 and study 2, respectively.
To evaluate the null distributions of the proposed metaanalysis statistics, we generated phenotype data sets under the null model. To reduce the computational burden, we simulated 20 replicates of covariates and phenotype for each of the 1000 genotype sets, which yielded 20,000 genotype-phenotype data sets. For scenario 1, the continuous phenotypes for the k-th (k¼ 1; 2; 3) study were generated by means of the following linear model: where sex k ¼ ðsex k1 ; Á Á Á ;sex kn k Þ T is a binary covariate vector in which each element sex ki is drawn from a Bernoulli distribution with probability 0.5; age k ¼ ðage k1 ; Á Á Á ;age kn k Þ T and bmi k ¼ ðbmi k1 ; Á Á Á ;bmi kn k Þ T are continuous covariate vectors in which each element is normally distributed, age ki $ Nð50;5 2 Þ and bmi ki $ Nð25;4 2 Þ, for 1˚i˚n k ; and e k ¼ ðe k1 ; Á Á Á ;e kn k Þ T represents the random errors, with each element following a standard normal distribution. For scenarios 2 and 3, we generated the phenotypes of study 1 according to Equation (18) by removing the age and gender effects, and generated the phenotypes of study 2 by removing the gender effect. The phenotypes of study 3 were generated with the full covariates. Under the null hypothesis, the genotypes are not associated with the phenotype.
To evaluate the statistical power of our proposed metaanalysis methods, we generated phenotypes under the alternative model based on the 1000 genotype data sets described previously. For each of the three scenarios, we considered five different levels of genetic main and interaction effects: level 1 refers to the existence of only genetic main effects but without interaction effects, level 2 refers to the existence of genetic main effects and weak interaction effects, level 3 refers to the existence of both genetic main effects and interaction effects, level 4 refers to the existence of interaction effects and weak genetic main effects, level 5 refers to the existence of only interaction effects but without genetic main effects. For each of the five levels, we assumed that 20, 40, or 60% of the rare variants were causal variants. We simulated one replicate of phenotype and covariates for each genotype data set. The covariates followed the same distributions described previously. For scenario 1, the phenotypes for the k-th study were generated by means of the following linear model: where G k is the genotype matrix of the causal variants in study k and b k has elements b kj $ Nð0;d 2 1 Þ. We used body mass index (BMI) as the environmental factor E k ¼ diagE ki , which was centered to have a mean of 0. c k represents the effects of gene-BMI interaction, with elements c kj $ Nð0;s 2 1 Þ. For scenarios 2 and 3, partial covariates were used as described in Table 1. The values of d 1 and s 1 for the five levels of genetic main and interaction effects were: level 1 with d 1 ¼ 0:2; s 1 ¼ 0, level 2 with d 1 ¼ 0:2; s 1 ¼ 0:005, level 3 with d 1 ¼ 0:2; s 1 ¼ 0:05, level 4 with d 1 ¼ 0:002; s 1 ¼ 0:05, level 5 with d 1 ¼ 0; s 1 ¼ 0:05.
For the case of homogeneous variant and interaction effects across studies, under all the three scenarios, we simulated the same variant effects and gene-BMI interaction effects across all studies, that is, b 1 ¼b 2 ¼ Á Á Á ¼b K ¼b and c 1 ¼c 2 ¼ Á Á Á ¼c K ¼c. For the case of heterogeneous variant and interaction effects across studies, under scenarios 1 and 2, we simulated the variant effects and gene-BMI interaction effects for each study independently. Under scenario 3, we simulated the same variant effects and gene-BMI interaction effects for study 1 and study 2 that have the EUR ancestry, and simulated the variant effects and gene-BMI interaction effects for study 3 with AA ancestry separately.
In all simulations, the variant weights followed Beta distributions, w j $ BetaðMAF j ; 1; 25Þ, as suggested by Wu et al. The gene-BMI interaction was considered significant if its Pvalue was less than the significance level of 2:5 Â 10 À6 , and empirical power was calculated as the proportion of significant results among 1000 replicates.

Null distributions of the meta-statistics
We examined the distributions of our proposed HOM-INT-FIX, HOM-INT-RAN, HET-INT-FIX, and HET-INT-RAN meta-statistics using the data sets simulated under the null hypothesis. We first computed the G Â E score statistics of each study according to formulas (5) and (7) and the regional G Â E relationship matrix of each study according to formulas (8) and (9). Then, for scenarios 1 and 2, we combined the study-level G Â E score statistics using formulas (10)- (11) and (14)-(15) to obtain the four meta-statistics. For scenario 3, we combined the study-level G Â E score statistics using formulas (10)- (11) and (16)-(17). We computed the regional G Â E relationship matrices for the meta-analysis by combining the regional relationship matrices from each study according to Equations (12) and (13). The P-values of the four meta-statistics were calculated according to their theoretical distributions and Kuonen's saddlepoint method. The distributions of the empirical P-values under the null hypothesis were compared with the expected uniform distribution between 0 and 1.
Quantile-quantile (Q-Q) plots of the meta-analysis statistics for testing G Â E effects under the three scenarios are shown in Figures 1-3, respectively. It can be observed that under all of the three scenarios, the empirical distributions of our four proposed statistics match well with their expected theoretical distributions.

Statistical power of the meta-statistics
The results for the statistical power under three scenarios are shown in Figures 4-6, respectively, which include results from five levels of main and interaction effects and three proportions of causal variants. In each figure, the power values on the left are computed based on the simulated data sets in which the variant and gene-BMI interaction effects were the same across studies. The results on the right are based on the data sets with heterogeneous variant and gene-BMI interaction effects across studies.
As can be seen that the statistical powers of the four proposed meta-statistics are close to 0 for interaction effects of levels 1 and 2 where there are no or weak interaction effects. For the interaction effects of levels 3-5, the powers are approximately the same while keeping meta statistics, simulation scenarios and proportions of causal variants fixed. Power results under scenario 1 are almost the same as scenario 2, which suggests that metaanalyses of studies with different covariates distributions had little influence on the power. This is because that the INT-FIX and INT-RAN statistics at study level are essentially based on the phenotypic residuals after adjusting the covariates. Statistical powers under scenario 3 are in general greater than those under scenario 1 or 2, which is due to that AA samples are included and average number of causal variants is larger than that under scenario 1 or 2.
It can be observed that HOM  Gene-age interactions in blood pressure traits with UK Biobank data UK Biobank is a prospective cohort study involving approximately 500,000 volunteers in the United Kingdom aged between 40 and 69 years, from whom extensive genetic and phenotypic data have been collected (Sudlow et al. 2015;Bycroft et al. 2018). The latest release of whole-exome sequencing data contains data from 200,643 participants. We extracted participants with European, African or Asian ancestry and excluded participants who had withdrawn and one member of each pair of first-or seconddegree relatives. The average SBP and diastolic blood pressure (DBP) measurements at the first visit were used. For participants taking antihypertensive medications at the time of the visit, 10 and 5 mm Hg were added to their SBP and DBP measurements, respectively (Cui et al. 2003). Mean arterial pressure (MAP) and pulse pressure (PP) were calculated as MAP ¼ SBP=3 þ 2DBP=3 and PP ¼ SBP À DBP, respectively, and the latter was then logarithmically transformed. Participants with missing BP phenotypes or covariates, including age, BMI and sex, were excluded. Phenotype and covariate outliers were also excluded, where these were defined as data points lying at least 5 standard deviations away from the corresponding means. In summary, a total of 162,148 European samples, 3113 African samples and 4745 Asian samples with complete phenotype, covariates, and exome genotype data were included in our analyses.
To conduct meta-analyses of multiple studies of the same ancestry, we divided the European samples into five groups based on the geographic locations of their assessment centers: group 1 included participants from assessment centers in Edinburgh, Glasgow, Middlesbrough, and Newcastle; group 2 was from Barts, Croydon, Hounslow, Oxford, and Reading; group 3 was from Bristol, Cardiff, and Swansea; group 4 was from Leeds, Sheffield, Nottingham, and Birmingham; and group 5 was from Bury, Cheadle, Liverpool, Manchester, Stockport, Stoke, and Wrexham.
For trans-ethnic meta-analyses, group 6 of African samples and group 7 of Asian samples were included. The characteristics of samples in the seven groups are presented in Table 2.
We selected nine loci that showed nominal evidence (P < 0.05) of SNP-age interaction in a genome-wide search of common variants with age-dependent effects in BP regulation (Simino et al. 2014). For the reported leading SNPs that are in gene regions, the corresponding genes were selected for analyzing gene-age interaction with rare variants; otherwise, the nearest up-and downstream genes were chosen. In total, 12 candidate genes were selected from the nine loci. The variants of the 12 genes were annotated with VEP (McLaren et al. 2016), and those annotated as stop_loss, missense_variant, start_lose, splice_    donor_variant, inframe_deletion, frameshift_variant, splice_ acceptor_variant, stop_gained, or inframe_insertion were used for analysis. In addition, variants that were PolyPhen or SIFT benign, defined as variants with PolyPhen scores smaller than 0.15 or SIFT scores larger than 0.05 (Cirulli et al. 2020), were excluded. PLINK (Chang et al. 2015) was used to extract rare variants with MAFs smaller than 0.03, and fcGENE (Roshyara and Scholz 2014) was used to convert genotypes into numeric values. The numbers of variants in the genes used for the single and multiple ancestries meta-analyses are shown in Tables 3 and 4   In Table 4

Discussion
In this study, we propose four meta-analysis methods for testing G Â E effects with rare variants while treating main genetic effects as either fixed or random and considering homogeneous or heterogeneous G Â E effects across studies. Simulations as well as an analysis of UK Biobank data demonstrate that treating variant main effects as either fixed or random provides approximately the same statistical power. The results are consistent with the conclusion in Chen et al. (2014) (Chen et al. 2014) at study level based on the analysis framework proposed by Lee et al. (2013). The methods offer much larger sample size for testing the interaction which is impossible for a single study. They were shown to be statistical efficient in our simulation studies as well as the analyses of UK Biobank data. It is well established that for the meta-analysis of association statistics with common variants, the power loss is minimal (Lin and Zeng 2010). However, the power loss of meta-analysis for rare variants is largely unexamined and has been suspected to be possibly more sizable (Panagiotou et al. 2013). In this study, we have compared the HOM-INT-FIX and HOM-INT-RAN meta-analyses with pooled analyses based on INT-FIX and INT-RAN, respectively. Our results show that they have approximately the same power. In gene-age interaction analyses of UK Biobank data, the results from the meta-analyses and the pooled analyses are close. Therefore, the power loss of meta-analysis for rare variants is concluded to be minimal as well.
In the gene-age interaction analysis of UK Biobank data, none of the 12 genes from the nine candidate loci shows experimentwide significant results. There are many possible reasons. First, most GWAS SNPs are from noncoding regions, and some of them may play regulatory roles by changing the expression levels of the modulated genes (Edwards et al. 2013). Nevertheless, the rare variants selected in our analyses likely alter the protein sequences that the genes express, which are not necessarily correlated with the regulatory SNPs. Second, it has been estimated that only approximately one-third of causal genes are the nearest genes to the GWAS loci (Gusev et al. 2016;Zhu et al. 2016). Thus, we may have missed the majority of the causal genes, which are located farther away from the candidate loci. Finally, the nine loci identified in the genome-wide analysis show only nominal evidence of interactions, which are not genome-wide significant. Therefore, some of the loci considered in the analyses may represent spurious interaction effects.
In this study, we have demonstrated the proposed meta-analyses for continuous traits. The methods can be readily generalized to the case of binary traits. For binary traits, logistic versions of INT-FIX and INT-RAN can be applied at the study level, and the G Â E score statistics for single variants and the relationship matrices can be computed based on (5) and (7)-(9), the same formulas as for continuous traits. There is no difference in how the meta-statistics are computed for binary and continuous traits. However, when studies have unbalanced case-control ratios and minor allele counts in a gene are very low, using saddlepoint approximation can result in inaccurate P-values for binary traits, and efficient resampling can be used instead .
Our proposed meta-analyses are limited to testing G Â E effects only. Joint testing of main genetic effects and interaction effects has long been suggested for the interaction analysis of common variants (Kraft et al. 2007;Manning et al. 2011). Joint testing offers better power than the analysis of main genetic effects only and the analysis of interaction effects only when both types of effects exist. In future work, we will further extend our metaanalyses to allow the joint testing of main genetic effects and G Â E effects.

Conclusions
In this study, we proposed four powerful meta-analysis methods for testing G Â E effects with rare variants. We considered both homogeneous G Â E effects and heterogeneous G Â E effects across studies and efficiently combined the summary statistics of INT-FIX and INT-RAN from multiple studies. Through simulations and real data analysis, we demonstrated that our approaches provide power comparable to that of pooled analysis when the interaction effects are homogeneous across studies. When heterogeneity exists across studies, our approaches can treat heterogeneity appropriately and achieve greater statistical power than a pooled analysis. Our meta-analysis methods of testing G Â E effects can be applied to synthesize results from multiple diverse studies to increase the effective sample size and improve the chance of identifying genes whose effects are modified by an environmental factor.

Funding
This work was supported by the national Thousand Youth Talents Plan.