{'fraction_non_alphanumeric': 0.044161135572057565, 'fraction_numerical': 0.02441348491818361, 'mean_word_length': 4.191914022517912, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 3, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SEAGLE: A Scalable Exact Algorithm for Large-Scale Set-Based GxE Tests in Biobank Data

The explosion of biobank data offers immediate opportunities for gene-environment (GxE) interaction studies of complex diseases because of the large sample sizes and the rich collection in genetic and non-genetic information. However, the extremely large sample size also introduces new computational challenges in GxE assessment, especially for set-based GxE variance component (VC) tests, which are a widely used strategy to boost overall GxE signals and to evaluate the joint GxE effect of multiple variants from a biologically meaningful unit (e.g., gene). In this work, we focus on continuous traits and present SEAGLE, a Scalable Exact AlGorithm for Large-scale set-based GxE tests, to permit GxE VC tests for biobank-scale data. SEAGLE employs modern matrix computations to achieve the same"exact"results as the original GxE VC tests without imposing additional assumptions or relying on approximations. SEAGLE can easily accommodate sample sizes in the order of $10^5$, is implementable on standard laptops, and does not require specialized computing equipment. We demonstrate SEAGLE's performance through extensive simulations. We illustrate its utility by conducting genome-wide gene-based GxE analysis on the Taiwan Biobank data to explore the interaction of gene and physical activity status on body mass index.


INTRODUCTION
Human complex diseases such as neurodegenerative diseases, psychiatric disorders, metabolic syndromes, and cancers are complex traits for which disease susceptibility, disease development, and treatment response are mediated by intricate genetic and environmental factors. Understanding the genetic etiology of these complex diseases requires collective consideration of potential genetic and environmental contributors. Studies of gene-environment interactions (G×E) enable understanding of the differences that environmental exposures may have on health outcomes in people with varying genotypes (Ottman, 1996;Hunter, 2005;McAllister et al., 2017). Examples include the impact of physical activity and alcohol consumption on the genetic risk for obesity-related traits (Sulc et al., 2020), the impact of air pollution on the genetic risk for cardio-metabolic and respiratory traits (Favé et al., 2018), and other examples reviewed in Ritz et al. (2017).
When assessing G×E effects, set-based tests are popular approaches to detecting interactions between an environmental factor and a set of single nucleotide polymorphism (SNPs) in a gene, sliding window, or functional region (Lin et al., 2013;Su et al., 2017;Lin et al., 2016;Tzeng et al., 2011;Zhao et al., 2015). Compared to single-SNP G×E tests, set-based G×E tests can enhance testing performance by reducing multiple-testing burden and by aggregating G×E signals over multiple SNPs that are of moderate effect sizes or of low frequencies.
Large-scale biobanks collect genetic and health information on hundreds of thousands of individuals. Their large sample sizes and rich data on non-genetic factors offer unprecedented opportunities for in-depth studies on G×E effects. While the explosion of biobank data collections provides great hopes for novel G×E discoveries, it also introduces computational challenges. In particular, many set-based G×E tests can be cast as variance component (VC) tests under a random effects modeling framework (Lin et al., 2013;Su et al., 2017), including kernel machine based tests (Lin et al., 2016;Wang et al., 2015a) and similarity regression based methods (Tzeng et al., 2011;Zhao et al., 2015). Hypothesis testing in this framework relies on computations with phenotypic variance matrices with dimension n × n (with n as the sample size) and may involve estimating nuisance variance components. When n is large, as in the case of biobank data, matrix computations whose operation counts scale with n 3 are prohibitive in terms of computation time and storage.
A number of methods attempt to ease this computational burden by bypassing the estimation of nuisance variance components, either through approximation of the variance or kernel matrices  or through approximation of the score-like test statistics (Wang et al., 2020). In the first case, approximating the kernel matrices still requires an expensive eigenvalue decomposition upfront, in addition to storage for the explicit formation of the n × n kernel matrices, thus lacking practical scalability. In the latter case, approximating the test statistics requires assumptions that may or may not be valid and are difficult to validate in practice. Our numerical studies in Section 3 show that the Type 1 error rates and power can be sub-optimal when data do not adhere to the required assumptions.
In this work, we focus on continuous traits and introduce a Scalable Exact AlGorithm for Large-scale set-based G×E tests (SEAGLE) for performing G×E VC tests on biobank data. Here "exact" refers to the fact that SEAGLE computes the original VC test statistic without any approximations, rather than the null distribution of the test statistic being asymptotic or exact. Exactness and scalability are achieved through the judicious use of modern matrix computations, allowing us to dispense with approximations and assumptions. Our numerical experiments illustrate that SEAGLE produces Type 1 error rates and power identical to those of the original GxE VC methods (Tzeng et al., 2011), but at a fraction of the speed. Additionally, SEAGLE can easily handle biobank-scale data with as many as n individuals in the order of 10 5 , often at the same speed as state-of-the-art approximate methods (Wang et al., 2020). Compared with the state-of-the-art approximate method (Wang et al., 2020), SEAGLE can produce more accurate Type 1 error rates and power. Another advantage of SEAGLE is its user-friendliness; it can be run on ordinary laptops and does not require specialized or high performance computing equipment or parallelization. In fact, our timing comparisons in Section 3 were performed on a 2013 Intel Core i5 laptop with a 2.70 GHz CPU and 16 GB RAM, specs that are standard for modern laptops. Therefore, SEAGLE makes it possible to run exact and scalable G×E VC tests on biobank-scale data with just a modicum of computational resources.
The rest of the paper proceeds as follows. Section 2 describes the standard mixed effects model for G×E effects, testing procedures, computational performance, and SEAGLE algorithm. Section 3 illustrates SEAGLE's performance through numerical studies. Section 4 concludes with a brief summary of our contributions and avenues for future work.

MATERIALS AND METHODS
We describe the standard mixed effects model for G×E effects and testing procedures (Section 2.1), the computational challenges for biobank-scale data (Section 2.2), the components of the SEAGLE algorithm (Section 2.3), and the SEAGLE algorithm as a whole (Section 2.3.4).

GxE Variance Component Tests for Continuous Traits
We present the standard mixed effects model for studying G×E effects, the score-like test statistic, and its p-value. Let y ∈ R n denote the response vector with n individual responses for a continuous trait; X ∈ R n×p the design matrix of p covariates whose leading column is the all ones vector for the intercept; E ∈ R n the design vector of the environmental factor in the G×E effect; and G ∈ R n×L the genetic marker matrix for the L SNPs where L < n. Define the design matrix for the G×E terms as G = diag(E)G ∈ R n×L where diag(E) ∈ R n×n is a diagonal matrix with the elements of the vector E on the diagonal.
Consider the linear mixed effects model (Tzeng et al., 2011;Lin et al., 2013), Here, β X ∈ R p and β E ∈ R are the fixed-effects coefficients for the covariates and environmental factor, respectively; b ∈ R L and c ∈ R L are the genetic main (G) effect and G×E effect, respectively, with b ∼ N(0, τ I L ) and c ∼ N(0, νI L ); ∼ N(0, σ I n ); and I k ∈ R k×k denotes the identity matrix of dimension k.
The SNP-set analysis models the G and G×E effects of the L SNPs as random effects rather than fixed effects. This choice avoids power loss for non-small L and numerical difficulties from correlated SNPs that can occur in a fixed effects model. To assess the presence of G×E effects with H 0 : c = 0 in Model (1), one can apply a score-like test to the corresponding variance component with H 0 : ν = 0.
To simplify the null model of (1) in the score-like test, we consolidate and define X = X E ∈ R n×P and β = β T X β T E T ∈ R P , where P = p + 1. The resulting null model becomes y = Xβ + Gb + , where the response is y ∼ N Xβ , V with V = τ GG T + σ I n . Following (Tzeng et al., 2011), the Frontiers score-like test statistic is In Equation (2) Appendix 1 presents the restricted maximum likelihood (REML) expectation-maximization (EM) algorithm for estimating the nuisance VC parameters τ and σ for computing T (Tzeng et al., 2011;Zhao et al., 2015). The test statistic T follows a weighted χ 2 (1) distribution asymptotically under H 0 : ν = 0. That is, T ∼ λ χ 2 (1) , where λ 's are the eigenvalues of Given the λ 's, the p-value of T can be computed with the moment matching method in Liu et al. (2009) or the exact method in Davies (1980).

Computational Challenges in GxE VC Tests for Biobank-Scale Data
We identify three computational bottlenecks.
1. The test statistic T and the p-value computation depend on P ∈ R n×n , which in turn depends on V −1 ∈ R n×n . Explicit formation of the inverse is too expensive and numerically inadvisable, due to loss of numerical accuracy and stability (Higham, 2002, Chapter 14).
2. The REML EM algorithm (Appendix 1) estimates the nuisance variance components τ and σ in V under the null hypothesis. Each iteration requires products with the orthogonal projector I − X( X T X) −1 X T , and inverting a matrix of dimension n − P ≈ n. 3. Computing the p-values requires two eigenvalue decompositions: (1) an eigenvalue decomposition of V to compute V 1 2 in C 1 ; and (2) an eigenvalue decomposition of C to compute the λ 's in the weighted χ 2 (1) distribution. Computing the eigenvalues and eigenvectors of the symmetric matrix V ∈ R n×n requires O(n 3 ) arithmetic operations and O(n 3 ) storage. Computing the eigenvalues of C ∈ R n×n requires another O(n 3 ) arithmetic operations.

Components of the SEAGLE Algorithm for Biobank-Scale GxE VC Test
We present our approach for overcoming the three computational challenges in the previous section: Multiplication with V −1 without explicit formation of the inverse (Section 2.3.1), a scalable REML EM algorithm (Section 2.3.2), and a scalable algorithm for computing the eigenvalues of C (Section 2.3.3). The idea is to replace explicit formation of inverses by low-rank updates and linear system solutions; and to replace n × n eigenvalue decompositions with L × L ones.

Multiplication by
The test statistic T and its p-value calculation depend on V −1 . We avoid the explicit formation of the inverse by viewing V = τ GG T + σ I n as the low-rank update of a diagonal matrix, and then applying the Sherman-Morrison-Woodbury formula below to reduce the dimension of the computed inverse from n to L where L n.
LEMMA 1 (Section 2.1.4 in Golub and Van Loan (2013)). Let H ∈ R n×n be nonsingular, and let U, B ∈ R n×L so that I + B T H −1 U is nonsingular. Then Applying Lemma 1 to the product of the inverse of V = σ I n + τ σ GG T with any right-hand side input W ∈ R n×l gives which reduces the dimension of the inverse from n to L. The explicit computation of the inverse of M = I L + τ σ G T G ∈ R L×L is, in turn, avoided with a Cholesky decomposition followed by a linear system solution. Algorithm 1 shows pseudocode for computing (4). As a further saving, we pre-compute the Cholesky factorization of M only once, so it is available for re-use in the computation of the test statistic and p-value.

Scalable REML EM Algorithm
We present the scalable version of the REML EM algorithm in Appendix 1 that avoids explicit formation of the orthogonal projector and the inverses.
We assume throughout that X has full column rank with rank( X) = P , and let range( X) be the space spanned by the columns of X. The space perpendicular to range( X) is range( X) ⊥ , and the orthogonal projector onto this space is where A ∈ R n×(n−P ) has orthonormal columns with A T A = I n−P . Let u = A T y ∈ R n−P be the orthogonal projection of the response onto range( X) ⊥ .
In iteration t + 1 of the algorithm, define

Frontiers
Then the updates in (S6) and (S5) from Appendix 1 can be expressed aŝ The two bottlenecks in the REML EM algorithm are the computation of the non-symmetric "square-root" A in (5), and products with R −1 from (6). Since P is small, n − P ≈ n, hence explicit formation of the inverse is out of the question, especially since R changes in each iteration due to the updates forτ t andσ t .
To avoid explicit formation of the full matrix in (5), we compute instead the QR decomposition where Q ∈ R n×n is an orthogonal matrix with Q T Q = QQ T = I n . The columns of Q 1 ∈ R n×P form an orthonormal basis for range( X), and the columns of Q 2 ∈ R n×(n−P ) form an orthonormal basis for range( X) ⊥ . The upper triangular matrix R 0 ∈ R P ×P is nonsingular, due to the assumption of X having full column rank. Therefore, (5) simplifies to Thus, A = Q 2 represents the trailing n − P columns of the orthogonal matrix Q in the QR factorization of X.
Algorithm 2 shows pseudocode for the REML EM algorithm. Since A = Q 2 occurs only as A T in matrix-vector or matrix-matrix multiplications, we do not compute A explicitly. Instead, we compute the full QR decomposition in (7) where the "QR object" Q = Q 1 A is stored implicitly in factored form.
To compute u = A T y, we multiplyũ = Q T y = Q T 1 y A T y and then extract the trailing n − P rows fromũ.
Furthermore, we apply R −1 from (6) with a modified version of Algorithm 1 where G is replaced by A T G and I n by I n−P . Unfortunately, one cannot pre-compute a Cholesky factorization for the whole algorithm sinceτ t andσ t change in each iteration. However, within a single iteration, we pre-compute a Cholesky factorization of R for subsequent linear system solutions of R −1 u and R −1 A T G. Following previous work on G×E VC tests in Tzeng et al. (2011), our convergence criteria are: (i) the magnitude of the relative difference between the current and previous estimate; and (ii) the default convergence tolerance from the SIMreg package for R.

Scalable Algorithm for Computing the Eigenvalues of C
Computation of the p-values requires the eigenvalues of C = C 1 C T 1 in (3), which in turn involves products with V 1 2 ∈ R n×n . We avoid the computation of the square root by exploiting the fact that the nonzero eigenvalues of C 1 C T 1 are equal to the nonzero eigenvalues of C T 1 C 1 . The symmetry of V and P and the equality PVP = P imply the much simpler expression The explicit formation of P is avoided by computing instead products P G with Algorithm 1. Therefore, our approach of replacing the n × n matrix C 1 C T 1 with the much smaller L × L matrix C T 1 C 1 reduces the operation count from O(n 3 ) down to O(L 3 ). Part III of Algorithm 3 shows the pseudocode.

The SEAGLE Algorithm
Combining the algorithms from Sections 2.3.1-2.3.3 gives the SEAGLE Algorithm 3 for computing the score-like test statistic T and its p-value. SEAGLE is implemented in the publicly available R package SEAGLE.
Algorithm 3 consists of three parts. Part I computesτ andσ with the scalable REML EM in Algorithm 2; Part II computes the score-like T statistic in (2); and Part III computes the p-values from the eigenvalues of C in (3). Linear systems with V are efficiently solved with Algorithm 1. The fast diagonal multiplication in R stores diagonal matrices as vectors. The QR decomposition is implemented with the qr function in the R base package. The qr.qty function makes it possible to left multiply by Q T without having to explicitly form Q.

Simulation Study
We evaluate the performance of our proposed method SEAGLE using simulation studies from two settings: I) data simulated from a random effects genetic model with n = 5,000 observations, and II) data simulated from a fixed effects genetic model with n = 20,000 and n = 100,000 observations. In the random effects simulations, we generate data according to Model (1). This enables us to evaluate SEAGLE's estimation and testing performance. We consider a smaller n to enable comparisons with existing G×E VC tests. In the fixed effects simulations, we generate data from a fixed effects model. This enables us to evaluate the testing performance when the data do not follow our modeling assumptions. We consider larger n values to demonstrate SEAGLE's effectiveness on biobank-scale data.
In each setting, we study the Type 1 error rate and power. We consider three baseline approaches: (i) the original G×E VC test (referred to as OVC) (Tzeng et al., 2011;Wang et al., 2015b), as implemented in the SIMreg R package (https://www4.stat.ncsu.edu/˜jytzeng/software simreg.php); (ii) fastKM , as implemented in the FastKM R package; and (iii) MAGEE (Wang et al., 2020), as implemented in the MAGEE R package. MAGEE is the state-of-the-art scalable G×E VC test with demonstrated superior performance compared to several set-based GxE methods.
In all simulations, we obtain the genotype design matrix G ∈ R n×L as follows. First, we employ the COSI software (Schaffner et al., 2005) to simulate 10,000 haplotypes of SNP sequences mimicking the European population. We then form a SNP set of L loci with minor allele frequency less than 1% by randomly selecting L SNPs without replacement. Finally, in each replicate, we generate the genotypes of n individuals by randomly selecting two haploytpes with replacement. We consider n = 5,000 and L = 100 in the random effects simulations, and n = 20,000 or 100,000 and L = 100 or 400 in the fixed effects simulations. We also consider a confounding factor X ∈ R n and an environmental factor E ∈ R n , where each is generated from a standard normal distribution. Given X and E, we then form the covariate design matrix X ∈ R n×3 by column-combining the vector of ones, X, and E together.

Random Effects Simulation Study
Given the genotype design matrix G of n = 5,000 and L = 100 loci and the n × P covariate design matrix X (where P = 3), we simulate the outcome data y according to the random effects model: y = Xβ + Gb + diag(E)Gc + e, where β is set as the all ones vector of length P ; b is generated from N(0, τ I L ); e is generated from N(0, σ I n ); and σ and τ are set to be 1. We set ν = 0 for Type I error analysis and ν > 0 for power analysis, where the actual value of ν is determined so that the empirical power is not too close to 1. We simulate N = 1,000 replicates and evaluate the results at the nominal level α = 0.05 for all analyses, except when assessing SEAGLE's Type I error rates at α = 5 × 10 −2 , 5 × 10 −3 and 5 × 10 −4 , where we consider N = 366,000.
We start by examining the Type 1 error rate for SEAGLE with N = 366,000 replicates. Table 1 shows that SEAGLE provides reasonable control over the Type 1 error rate at varying α-levels. Next, under H 0 : ν = 0 and τ = σ = 1 with N = 1,000 replicates, we compare the testing results of SEAGLE with OVC. Table 2 shows the bias and the mean square error (MSE) of the estimated values for τ and σ obtained from the SEAGLE and OVC REML EM algorithms. Both algorithms produce very small bias and MSE for τ and σ. Figures 1 (A) and (B) depict scatter plots of the score-like test statistics and p-values, respectively, produced by SEAGLE and OVC. The figures show that SEAGLE and OVC produce identical test statistics and p-values, hence the "exactness" of the SEAGLE algorithm. Since the data are generated from a random effects model under H 0 : ν = 0, we can compute the "true" score-like test statistic T by evaluating (2) at the true τ and σ values, and obtaining the corresponding p-value. We refer to this as the "Truth" and include it as a baseline approach. Supplementary Figure S1 depicts quantile-quantile plots (QQ plots) of the p-values obtained over N = 1,000 replicates for Truth, SEAGLE, OVC, FastKM, and MAGEE. All methods exhibit similar p-value behavior, and the red line for SEAGLE is very close to the light blue (OVC) and yellow lines (Truth). In Table 3, we also compute the MSE of the p-values obtained from SEAGLE, OVC, FastKM, and MAGEE, compared to the Truth p-values. We observe that MAGEE produces p-values with larger MSE than the other methods. Supplementary Figure S2 shows the corresponding absolute relative error of the p-values for each method, computed by first taking the absolute difference between a method's p-value and the Truth p-value, then dividing it by the Truth p-value. The boxplots suggest that MAGEE exhibits higher bias and greater variance than SEAGLE, OVC and FastKM.
Regarding the computational cost, Figure 2 (A) shows boxplots of the computation time in seconds required to obtain a single p-value for each of the methods over N = 1,000 replicates with τ = σ = 1 and ν = 0. Figure 2 Figure 3 shows the Type 1 error rate for each method under H 0 : ν = 0. SEAGLE performs identically to OVC with respect to Type 1 error rate at α = 0.05 while requiring only a fraction of the computation time, as demonstrated in Figure 2. By contrast, MAGEE is nearly as fast as SEAGLE for n = 5,000 and L = 100 but produces much more conservative p-values at τ = σ = 1. Figure 4 shows the power for each method under the alternative hypothesis that ν > 0. SEAGLE again performs identically to OVC while requiring a fraction of the computation time. By contrast, MAGEE is nearly as fast as SEAGLE for n = 5,000 and L = 100 but has lower power when τ = σ = 1 and ν = 0.04.

Fixed Effects Simulation Study
To study the performance of our proposed method when the data may not adhere to our model assumptions, we follow previous work Wang et al., 2020) and simulate data according to the fixed effects model with a given X and G: y = Xγ X + Gγ G + diag(E)Gγ GE + e, where γ X is the all ones vector of length P = 3, γ G ∈ R L , γ GE ∈ R L , and e ∼ N(0, σ I n ). The entries of γ G and γ GE pertaining to causal loci are set to be γ G and γ GE , respectively. The remaining entries of γ G and γ GE pertaining to non-causal loci are 0. We consider n = 20,000 or 100,000 observations with L = 100 or 400. We select the first loci to be causal (i.e., loci with both G and G×E effects or just G effect). We vary γ G over 0.5, 1, and 1.5 for the loci to study the impact of the G main effect sizes. We compare SEAGLE with MAGEE only since OVC and fastKM are unable to work on the sample sizes considered here.
We first evaluate the Type 1 error of SEAGLE by simulating N = 1,000 replicates with L = 100 for both n = 20,000 and 100,000, and setting γ GE = 0 for all loci while letting the first = 40 loci to have non-zero γ G . Figure 5 depicts the Type 1 error rate at α = 0.05 over varying values for γ G . While the Type 1 error rate for SEAGLE remains relatively unaffected by different γ G values, MAGEE produces more conservative p-values as γ G increases. This is consistent with the MAGEE assumption requiring a small G main effect (Wang et al., 2020). Supplementary Figure S3 shows the corresponding quantile-quantile plots for the p-values obtained from SEAGLE and MAGEE. Figure 6 depicts boxplots of the computation time in seconds required to obtain a single p-value over the N = 1,000 replicates for n = 20,000 and n = 100,000, and L = 100, over varying values for γ G . All replicates were computed on a 2013 Intel Core i5 laptop with a 2.70 GHz CPU and 16 GB RAM. For n = 20,000, SEAGLE is faster than MAGEE at larger values of γ G even though MAGEE computes an approximation to the test statistic T and bypasses the traditional REML EM algorithm. At smaller values of γ G , however, SEAGLE requires a few seconds more than MAGEE. This is because smaller γ G values  are determined so that the power for n = 20,000 at α = 0.05 is not close to 1. The values of γ G for the causal loci are set to be 0.5, 1.0, and 1.5 as before. Figure 7 shows the power for L = 100 loci. At n = 20,000, SEAGLE exhibits better power than MAGEE at all combinations of γ G and γ GE . Moreover, the difference in power increases for larger values of γ G since MAGEE relies on the assumption that the G main effect size is small. At n = 100,000 and the same values of γ GE , we report the power at α = 0.001 instead of 0.05 because the power at α = 0.05 is near 1 for both methods. We see that both methods produce similar results although SEAGLE still outperforms MAGEE at slightly smaller values of γ GE . Figure 8 shows the power for L = 400 loci. Similar patterns of relative power performance are observed as in the case of L = 100, except that the power difference between SEAGLE and MAGEE is more pronounced in L = 400.

Application to the Taiwan Biobank Data
To illustrate the scalability of the G×E VC test using SEAGLE, we apply SEAGLE and MAGEE to the Taiwan Biobank (TWB) data. TWB is a nationwide biobank project initiated in 2012 and has recruited more than 15,995 individuals. Peripheral blood specimens were extracted and genotyped using the Affymetrix Genomewide Axiom TWB array, which was designed specifically for a Taiwanese population. We conduct the gene-based G×E analysis and evaluate the interaction between gene and physical activity (PA) status on body mass index (BMI), adjusting for age, sex and the top 10 principal components for population substructure. The PA status is a binary indicator for with/without regular physical activity. Our G×E analyses focuses on a subset of 11,664 unrelated individuals who have nonmissing phenotype and covariate information. After PLINK quality control (i.e., removing SNPs with call rates < 0.95 or Hardy-Weinberg Equilibrium p-value < 10 −6 ), there are 589,867 SNPs remaining, which are mapped to genes according to the gene range list "glist-hg19" from the PLINK Resources page at https://www.cog-genomics.org/plink/1.9/resources. There are a total of 13,260 genes containing >1 SNPs for G×PA interaction analysis.
The median run time of SEAGLE and MAGEE is 2.4 and 1.3 seconds, respectively, Both SEAGLE and MAGEE do not find any significant G×PA interactions at the genome-wide Bonferroni threshold 0.05/13,260 = 3.77 × 10 −6 . We hence discuss the results using a less stringent threshold, i.e., 5 × 10 −4 and summarize the results in Figure 9 and Supplementary Table S1. SEAGLE and MAGEE identify 8 and 6 G×PA interactions, respectively, among which 5 G×PA results are identified by both methods (Figure 9). The observation that SEAGLE identifies slightly more G×PA effects than MAGEE generally agrees with the simulation findings. We use the GeneCards Human Gene Database (www.genecards.org) (Stelzer et al., 2016) to explore the relevance of the identified genes with BMI or PA (see Supplementary Tables S1). Two of the 5 commonly identified genes, i.e., FCN2 and OCM, have non-zero relevance scores (i.e., 0.56 and 0.91, respectively). For the 3 genes identified by SEAGLE only, i.e., ALOX5AP, BCLAF1 and PCDH17, their relevance scores are 6.16, 0.26 and 1.54, respectively. The expression of ALOX5AP has also been found associated with obesity and insulin resistance (Kaaman et al., 2006) as well as exercise-induced stress (Hilberg et al., 2005). On the other hand, TBPL1 (identified by MAGEE only) is not in the GeneCards relevance list with BMI or PA.

DISCUSSION
We introduced SEAGLE, a scalable exact algorithm for performing set-based G×E VC tests on large-scale biobank data. We achieve scalability and accuracy by applying modern numerical analysis techniques, which include avoiding the explicit formation of products and inverses of large matrices. Our numerical experiments illustrate that SEAGLE produces Type 1 error rates and power that are identical to those of the original VC test (Tzeng et al., 2011), while requiring a fraction of the computational cost. Moreover, SEAGLE is well-equipped to handle the very large dimensions required for analysis of large-scale biobank data.
State-of-the-art computational approaches such as MAGEE bypass the traditional time-consuming REML EM algorithm, and instead compute an approximation to the score-like test statistic by assuming that the G main effect size is small. In practice, however, the G main effect size is often unknown. Our numerical experiments illustrate that SEAGLE generally achieves better Type 1 error and power with comparable computation time.
We highlight the fact that our timing experiments were performed on a 2013 Intel Core i5 laptop with a 2.70 GHz CPU and 16 GB RAM. Therefore, SEAGLE performs efficient and exact set-based G×E tests on biobank-scale data with n = 20,000 and n = 100,000 observations on ordinary laptops, without any need for high performance computational platforms. This makes SEAGLE broadly accessible to all researchers. Software for SEAGLE is publicly available as the SEAGLE package on GitHub (https://github.com/jocelynchi/SEAGLE), and will soon be available on the Comprehensive R Archive Network as well.
We conclude with a discussion of three avenues for future extensions. First is the extension of Model (1) to generalized linear mixed models that can accommodate binary phenotype traits in addition to the continuous phenotype covered here.
Second is the extension from a single environmental factor to a set of factors represented by E ∈ R n×q with q > 1. The corresponding extension of Model (1) is h E ∼ N (0, τ E Σ E ), h G ∼ N (0, τ G Σ G ), h GE ∼ N (0, νΣ GE ), and ∼ N (0, σI n ).
Here Σ E , Σ G , Σ GE ∈ R n×n are variance matrices, where Σ E = EE T , Σ G = GG T and Σ GE is the element-wise products of Σ G and Σ E . A straightforward adaptation of SEAGLE's scalable REML EM algorithm to the EM algorithm in Wang et al. (2015b) takes care of estimating the nuisance VC parameters. Numerical analysis techniques analogous to the ones presented here will be the foundation for the efficient extension to multi-E factors. Table 1. Type 1 error rate with 95% confidence intervals (CIs) for SEAGLE over N = 366,000 replicates for n = 5,000 observations, L = 100 loci, and variance components τ = σ = 1 and ν = 0.