Subset-Based Analysis Using Gene-Environment Interactions for Discovery of Genetic Associations across Multiple Studies or Phenotypes

Abstract

Objectives: Classical methods for combining summary data from genome-wide association studies only use marginal genetic effects, and power can be compromised in the presence of heterogeneity. We aim to enhance the discovery of novel associated loci in the presence of heterogeneity of genetic effects in subgroups defined by an environmental factor. Methods: We present a pvalue-assisted subset testing for associations (pASTA) framework that generalizes the previously proposed association analysis based on subsets (ASSET) method by incorporating gene-environment (G-E) interactions into the testing procedure. We conduct simulation studies and provide two data examples. Results: Simulation studies show that our proposal is more powerful than methods based on marginal associations in the presence of G-E interactions and maintains comparable power even in their absence. Both data examples demonstrate that our method can increase power to detect overall genetic associations and identify novel studies/phenotypes that contribute to the association. Conclusions: Our proposed method can be a useful screening tool to identify candidate single nucleotide polymorphisms that are potentially associated with the trait(s) of interest for further validation. It also allows researchers to determine the most probable subset of traits that exhibit genetic associations in addition to the enhancement of power.

© 2019 S. Karger AG, Basel

Introduction

Genome-wide association studies (GWAS) are popular epidemiologic tools for studying the genetic architecture underlying a phenotypic trait [1]. Meta-analysis is a commonly used approach to combine genetic associations across multiple independent studies [2]. Fixed-effect meta-analysis based on summary statistics is known to retain full efficiency of an analysis based on individual level data [3]. Since fixed-effect meta-analysis assumes common effect sizes across studies, the power of the fixed-effect method may be compromised when there is between-study heterogeneity in the effect size of the association. Han and Eskin [4] have proposed a random-effect model with a likelihood ratio test to account for such heterogeneity. More recently, methods for aggregating association signals across multiple related phenotypes have also become more popular [5, 6].

The association analysis based on subsets (ASSET) test, proposed by Bhattacharjee et al. [7], uses summary statistics from individual association analyses to develop a powerful test that allows the existence of both a subset of null results and effects in opposite directions in different individual tests. Unlike the random-effect model proposed by Han and Eskin [4], ASSET assumes fixed effects for individual associations but still allows for heterogeneity in the data. This approach essentially explores all possible non-empty subsets of available studies/traits and searches for the one that yields the strongest evidence of association while adjusting for multiple testing during such exhaustive search. One appealing feature of ASSET is that the most probable subset of studies/traits that exhibits genetic association can be identified in addition to the enhancement of power of detecting association signals. This method has been widely used in analyses of associations across multiple studies and phenotypes [8, 9].

Exploiting gene-environment (G-E) interaction may help identify genetic variants that do not demonstrate very strong marginal effects due to environmental heterogeneity but may have stronger genetic effects in subgroups defined by certain levels of exposure [10-13]. Like most commonly used methods that screen only for marginal associations, ASSET does not account for potential G-E interactions in the testing procedure. If a genetic variant only affects a subgroup, for example the exposed group, then in the presence of such pure G-E interaction, the power of ASSET can be largely compromised. Several existing versions of 2-degree-of-freedom (2-df) tests for discovering genetic associations take G-E interactions into account. Kraft et al. [12] proposed a likelihood ratio test that jointly tests for genetic main effects and G-E interactions in case-control studies. Dai et al. [11] proposed a slightly different approach that is based on the sum of two Wald χ² statistics for detecting marginal effects combined with G-E interactions. The latter approach has the advantage of incorporating G-E independence into the test statistic for G-E interaction in case-control studies to achieve increased power [14, 15]. The literature on meta-analysis of G-E interactions or joint tests remains relatively limited [16-18], and none of these approaches attempt to specifically identify a subset of studies/traits that is most likely to have genetic associations and/or interactions. Current literature shows that the 2-df tests mentioned above can offer enhanced power in the presence of G-E interactions and do not lead to significant loss of power even when the interactions are absent [16, 18].

In this paper, we propose a modification to ASSET that incorporates G-E interactions to increase the power of detecting genetic associations in combining summary data across studies of a given phenotype (say type 2 diabetes, T2D) or across multiple related phenotypes (say a set of lipid traits) measured in different studies. Our approach is similar to ASSET in terms of searching over all possible combinations of studies/phenotypes and, thus, inherits its advantage of identifying the maximally associated subset along with increased power for detecting the overall associations. The proposed framework uses p values from the underlying tests as input instead of Z-statistics that are used by ASSET, and it is, therefore, referred to as pvalue-assisted subset testing for associations (pASTA). To be able to incorporate G-E interactions, pASTA uses 2-df tests as the underlying test statistics to be combined. In fact, the p values can result from any statistical tests, including tests with multiple degrees of freedom. On the other hand, by using p values, one loses the ability to incorporate the signs of coefficients (for both genetic main effect and G-E interaction), and the method is intrinsically two-sided. Since the key input variables are the p values from underlying tests, pASTA can be used under any commonly used epidemiologic study design, and we present our simulation studies and data examples for both case-control and cohort studies. Like ASSET, the analytical method in pASTA can handle a set of correlated p values, enabling its validity in studies with overlapping subjects or multiple phenotypes measured on the same set of subjects.

The rest of the paper is organized as follows. We first describe the construction of our test statistic and derivation of the associated p value. We present extensive simulation studies to compare the performance of our proposed method to ASSET and Fisher’s combined p value approaches [19] with and without G-E interactions under various parameter settings. Specifically, we examine the type I error rates and the power to detect (1) the truly associated variant and (2) the truly associated subset of studies/traits for the given variant. We illustrate our proposed method with two data applications. The first is a meta-analysis of 6 case-control studies of T2D among European populations (3 are Finnish, 1 is German, and 2 are Norwegian) with 4,422 cases and 5,202 controls. We focus on two single nucleotide polymorphisms (SNPs) in the FTOgene (MIM 610966) related to obesity and their interactions with body mass index (BMI). The second application is a multiple-phenotype analysis of 9 lipid-related quantitative traits using data of 5,123 individuals from the North Finland Birth Cohort 1966 (NFBC1966) study. We investigate the association between two SNPs discovered by GWAS, one near the LPL gene (MIM 609708) and the other in the APOBgene (MIM 107730), and the 9 lipid traits while taking the SNP-BMI interactions into consideration. We conclude the paper with a discussion.

Material and Methods

Notation

We consider K studies, each with sample size n_k. We allow the studies to have a set of overlapping subjects. Let Y_ki, G_ki, and E_ki denote the underlying phenotype, the given genetic variant, and the environmental factor, respectively, measured for the i-th subject in the k-th study. The trait(s) Y can either be disease status in a case-control study or binary or quantitative trait(s) measured in a cohort study. The genetic variant G may be coded as binary (under a dominant or recessive susceptibility model) or as allele count (under the additive model), and for the latter case, we treat G as continuous dosage. The environmental factor E can either be categorical or continuous. We mainly consider two models for each of the K studies: the marginal genetic association model

and the joint model with G-E interaction

where i = 1, …, n_k and k = 1, …, K. The choice of the link function g(•) depends upon the type of Y. Adjusted covariates are dropped from the presentation for simplicity of notations. Note that the K studies in this setup can be replaced by K traits in a single study of size n, and the inferential procedure for each of the K traits stays the same with n_k = n for each k = 1, …, K.

We primarily consider three types of association tests for each of the K studies: (1) marginal genetic effect in Model 1 (MA), with null hypothesis (H_0k) being α^(k)_G = 0; (2) genetic main effect and G-E interaction in Model 2 (JOINT) [12], with H_0k being β^(k)_G = β^(k)_GE = 0; and (3) marginal genetic effect in Model 1 and G-E interaction in Model 2 (MA+GE) [11], with H_0k being α^(k)_G = β^(k)_GE = 0. MA is a 1-df test, while JOINT and MA+GE are 2-df tests. A summary description of the different tests and their corresponding null hypotheses, test statistics, and distributions under the null hypothesis are presented in Appendix Table A1. The test of the form MA+GE is typically the sum of two Wald χ² test statistics that are known to be asymptotically independent and can, thus, be combined to yield a 2-df χ² statistic [20].

For a case-control study, the coefficient of interaction β^(k)_GE can be estimated using case-control (CC), case-only (CO) [14], or empirical Bayes (EB) [15] estimators, denoted as β^^(k)_CC, β^^(k)_CO, and β^^(k)_EB, respectively. The CO estimator increases the power of detecting G-E interaction by using G-E independence but incurs inflated type I error under violations of such independence. The EB estimator can be viewed as a linear combination of CC and CO estimators, and it improves power to detect G-E interaction when compared to standard CC analysis and provides better control of type I error when compared to a CO analysis. It has been proved that β^^(k)_CC and α^^(k)_G are asymptotically independent, and so are β^^(k)_EB and α^^(k)_G when we assume G-E independence and a rare disease [20]. As such, for case-control studies, we use β^^(k)_CCand β^^(k)_EB separately as two possible choices for the estimator of β^(k)_GE when computing the test statistic of MA+GE and denote these two types of tests as MA+CC and MA+EB, respectively [11]. Note that the use of G-E independence requires the use of a retrospective likelihood of the form (G, E)│D, which is not valid unless the sampling is conditional on D. A cohort analysis is always based on D│(G, E), and we cannot model the stochastic distribution of (G, E). Therefore, for a cohort study with multiple quantitative or binary phenotypes, the use of G-E independence does not lead to enhanced power and as such the EB estimator does not apply.

p Value-Assisted Subset Testing for Associations

Regardless of the test statistic chosen for each study, the underlying analytical framework of pASTA starts with a set of p values: p_k, k = 1, …, K. In our setting, they are obtained from the three tests MA, JOINT, or MA+GE. For pASTA, we define the Z-statistic for the k-th study as

Z_k = –Φ^–1(p_k),

where Φ(•) is the cumulative distribution function corresponding to the standard normal distribution. That is, smaller p values indicate larger Z-statistics, and the null hypothesis H_0k should be rejected if Z_k is large in the positive direction. For any non-empty subset S ⊆ {1, …, K}, let

where we consider

as the weight of Z_k. Other options for the weights

can also be accommodated in this framework. For example, when cases and controls are not of equal numbers, n_k can be substituted by the effective sample size [2]. The test statistic for evaluating the overall association of a given genetic variant is defined as

and is obtained by searching over all possible │𝓢│ = 2^k – 1 non-empty subsets of studies.

It is challenging to characterize the exact distribution of Z_meta in a closed form, especially when the number of studies, K, is large. Therefore, we borrow the idea followed in ASSET and use the discrete local maxima (DLM) procedure [21] to obtain a conservative upper bound for the meta-analytic combined p value. The DLM procedure automatically accounts for the multiple comparison issue. The detailed derivation of the combined p value based on Z_meta is included in the Appendix. In brief, let γ ∈ Γ be a K-dimensional vector, each of whose components takes value in {0, 1}, and let S_γ be the subset of studies whose corresponding coordinates in γ are 1. A neighbor of S_γ is defined as a subset in  that can be obtained by adding or dropping one study to or from S_γ. Then, Z(S_γ^*) is called a local maximum, if it is greater than all Z(S_γ)s, where S_γ refers to any possible neighbor of S_γ^*. For an observed test statistic Z_meta = T_obs, the DLM-approximated p value can be expressed as (see also equation A1 in the Appendix)

where h(•) is the probability distribution function of Z(S_γ). This DLM-based p value is actually an upper bound to the exact one (equation A1 in the Appendix). Hence the overall procedure is conservative in terms of type I error.

Each of the terms of the form P(Z(S_γ ± k) < z│Z(S_γ) = z) can be evaluated based on the conditional normal distribution of Z_k given Z(S_γ) (equations A2 and A3 in the Appendix). Thus, all we need to evaluate is the joint bivariate distributions of Z_k and Z(S_γ) for k = 1, …, K. Note that Z_k follows the standard normal distribution under the null hypothesis that there is no association (and interaction) in study k. The set of Z_ks are independent if the studies are independent. For multiple studies with overlapping subjects or multiple phenotypes measured on the same set of individuals, the Z_ks will be correlated with variance 1 and correlation matrix Σ with entries, say, ρ_kk′ = Corr(Z_k, Z_k′).

Let L = (L₁, …, L_K), where L_i = 1 if the i-th study is in S_γ and L_i = 0 otherwise. Let Z = (Z₁, …, Z_K) follow a multivariate normal distribution MVN_K(0, Σ). In particular, if all studies/traits are independent, then Σ is a K × K identity matrix. Let w = (w₁, …, w_K) be the vector of weights for the element in L.

if L_i = 1 and 0 if L_i = 0. Let e_kbe a vector of length K where the k-th element equals 1 and 0 elsewhere, e.g., e₂ = (0, 1, 0, …, 0)^T. Let

which is a 2 × K matrix. Under the null hypothesis,

Several methods have been proposed to estimate the correlation matrix Σ. Lin and Sullivan [22] provided a formula to estimate the correlation between a pair of studies using the information of shared cases and controls when the number of overlapping subjects is known. They also discussed a modified version of the formula for quantitative traits in a cross-sectional study. Bhattacharjee et al. [7] extended Lin and Sullivan’s formula for case-control studies to accommodate the situation where cases of one study serve as controls of another. However, these estimates will likely be inaccurate when the overlap fraction is large, and therefore do not directly apply to multiple-phenotype analyses, where the phenotypes are measured on the same set of subjects and the overlap fraction is 100%. In this case, since a large positive correlation between phenotypic measurements usually indicates a large positive correlation between Z-statistics, an intuitive way to approximate the correlations of Z_ks is to use the phenotypic correlations [23], as is done in ASSET. We adopt this approach for our simulation studies and the second data example to estimate the correlation matrix Σ corresponding to Z. In other words, if

then we let Σ^ = S.

Comparison of pASTA and ASSET

We would like to point out a few features of pASTA as compared to ASSET when one is only interested in marginal genetic association estimated by Model 1. First, ASSET and pASTA use the Z-statistics in different ways (see Appendix Table A2). The Z-statistics used in pASTA no longer differentiate between possible directions of association, while the signs of the Z-statistics used in ASSET are consistent with those of α^^(k)_Gs. ASSET can detect genetic associations in either one direction (one-sided) or both directions (two-sided). For the latter case, ASSET searches for the subsets of studies that yield the strongest evidence separately in both directions. It obtains a p value for each direction, and it combines the two p values using Fisher’s combined p value method [19]. pASTA does not have this feature. A more detailed comparison of the Z-statistics, test statistics, and analytic expression for the DLM-based p values of the two methods are presented in Appendix Table A2. The extension in this paper is mainly to incorporate interaction while testing association using 2-df tests in the ASSET framework via the use of p values.

Simulation Studies

Simulation Study 1: Combining Signals across Multiple Independent Case-Control Studies

Methods Considered.We first assess the performance of the proposed pASTA approach for combining summary results from multiple independent case-control studies that have no overlapping subjects or relatedness. We compare pASTA to three classes of alternatives (Table 1): ASSET, Fisher’s combined p value method [19], and the gold standard test which assumes that the true subset of non-null studies and the true model from which the data are generated are known a priori, to reflect the maximum achievable power for benchmarking each method. Given K independent p values p₁, …, p_K, the χ² test statistic of Fisher’s method follows a χ²_{df= 2K} distribution under the global null. Bhattacharjee et al. [7] has already compared an inverse variance-weighted meta-analysis with ASSET and Fisher’s method, showing that the former is not as powerful as the latter two. Therefore, we do not include a standard fixed-effect inverse variance-weighted meta-analysis method in our comparative study unless otherwise noted. We use several choices of 2-df tests as described in Table 2 and expanded in Appendix Table A1.

Table 1.

Methods considered in simulation studies

Table 2.

Type I error rates at various significance levels for the 2 simulation studies

Data Generative Model. We consider K = 10 independent studies with 6,000 cases and 6,000 controls. Disease status D, genetic marker G, and environmental factor E are assumed binary (1 for presence and 0 for absence) for simplicity. Following the standard approach of simulating case-control data with G-E interaction [10], we generate the data from a logistic regression model

for k = 1, …, K. In addition, we let Ɵ_GE, p_G, and p_E denote the G-E association, the frequency of G = 1 (this probability depends on the frequency of the risk allele, we only consider a dominant model in the simulation study), and the prevalence of environmental factor (E = 1) in controls, respectively. We specify the parameters Ɵ_GE, p_G, p_E, β^(k)_G, β^(k)_E, and β^(k)_GE, and determine the odds ratios of marginal genetic association (OR_G) and environmental association (OR_E) using the formula in Boonstra et al. [10]. In all of our simulation studies, we assume G-E independence among the control population, i.e., Ɵ_GE =0.

Evaluation of Type I Error. To assess the type I error rate of each method in the presence or absence of environmental effects, we generate data under the null hypothesis of no genetic effects, namely, set β^(k)_G= β^(k)_GE = 0 in Model 2. We consider two scenarios (Appendix Table A3), one with no environmental factor and hence no G-E interaction and the other with OR_E =1.4. All ten studies share the same parameters in each scenario (Appendix Table A3). We evaluate the type I errors at three levels of significance: α = 10^–2 and 10^–3 resembling a candidate gene study, and α = 10^–5 resembling a large-scale exploration. We simulate 5 × 10⁵ replicated datasets for each scenario.

Metrics for Power Analysis. We compare the power of pASTA with the other methods displayed in Table 1 from two perspectives: the power of detecting overall genetic associations and the accuracy of identifying the exact subset of non-null studies. Here, the non-null studies refer to studies in which either α^(k)_G ≠ 0 (and hence β^(k)_G ≠ 0) or β^(k)_GE ≠ 0. The specific definitions of the two metrics are given in the Appendix. We also report the sensitivity and specificity of identifying the subset of non-null studies, defined as the proportion of non-null studies that are correctly identified by subset-based approaches, and the proportion of truly null studies that are declared null, respectively.

Power for Detecting Marginal Association. Prior to consideration of G-E interaction, we extensively investigate the relative performance of ASSET and pASTA(MA) in detecting marginal genetic associations in a case-control setting where the data are generated under Model 1. We consider two scenarios in this simulation, each with K = 10 independent studies and k₀ = 4 of them having true associations. The prevalence of the genetic variant p_G is 0.2. In the first scenario, we assume that the marginal genetic effects in the 4 non-null studies are all in the positive direction. In the second, we reverse the sign of the log odds ratio parameter of genetic effect (β_G) in 2 of the 4 studies. In both scenarios, the absolute values of the log odds ratio of marginal associations │α^(k)_G │ are the same for all non-null studies, ranging from log 1.05 to log 1.2, and OR^(k)_G is incremented by 0.01. Since in reality it is possible that all studies have non-null signals, whether constant or varying across studies, we also carry out a simulation study where k₀ = 10. For this particular simulation setting, we include a fixed-effect meta-analysis method in addition to the 1-df tests described in Table 1. For all the settings described above, we use α = 0.001 as the significance level, which resembles a candidate gene study with 50 SNPs, and simulate 2,000 datasets for each parameter configuration.

Power for Detecting Association after Incorporating Interaction.We then consider both G and E in our model and the whole set of 1-df and 2-df tests as listed in Table 1. We consider 5 scenarios (Appendix Table A4) to evaluate the power of pASTA. Scenario 1 assumes no G-E interaction and a moderate marginal genetic effect (OR_G). In Scenario 2, we set parameters to induce large G-E interaction but small marginal genetic effect. Scenario 3 considers the same magnitude of interaction as Scenario 2, but with zero genetic and environmental main effects, which result in a smaller marginal genetic effect compared to Scenario 2. Scenario 4 assumes strong G-E interaction and protective genetic main effect, the two of which cancel out and result in an odds ratio of marginal genetic effect close to 1. Scenario 5 considers negative interaction and positive genetic main effect, which lead to a nearly-null marginal genetic effect. The last 2 scenarios represent less common case-control settings, where the genetic associations exist in opposite directions for exposed and unexposed groups. The number of non-null studies ranges from 2 to 7. Simulations are repeated 2,000 times for each scenario, and the significance level is set at α = 0.001.

Simulation Study 2: A Multiple-Phenotype Analysis in a Cohort Study

Methods Considered.We compare the performance of pASTA and ASSET in a cohort study with a set of multiple correlated traits. In this case, the test statistics of Fisher’s method no longer follows a χ²_df=2K distribution under the null hypotheses due to the correlations among p₁, …, p_K. Resampling-based approaches can be used to obtain the null distribution when correlations exist [24], but in this situation, we are dealing with a significance level at the magnitude of 10^–5 or smaller. It is computationally intensive to obtain a null distribution to evaluate the p value at such a stringent threshold. Therefore, we exclude comparison with Fisher’s method from Simulation Study 2.

Data Generative Model.We consider a multiple-phenotype study with 6,000 subjects and K = 9 correlated traits. The outcome Y is generated from a multivariate linear regression model

where (ε_1i, …, ε_Ki)^T ∼ MVN(0, Ω). We use the empirical phenotypic correlations from the NFBC1966 study (Appendix Table A5) as the generative covariance matrix Ω. We also assume a compound symmetry structure for the correlation matrix and vary the correlation (ρ = 0.2, 0.5, or 0.7) for the power analysis.

Type I Error and Power Analysis.Similar to Simulation Study 1, the type I error rates are estimated in the presence and absence of an environmental factor (Appendix Table A6). The power of these subset-based approaches are evaluated in 2 scenarios (Appendix Table A7) using the same metrics as used in Simulation Study 1. In Scenario 1, no interaction is involved, while in Scenario 2, interaction comes into play.

Data Analyses

We illustrate our methods by using data from 2 genetic association studies. The first is a meta-analysis of 6 independent case-control studies of T2D and the other considers a multiple-phenotype analysis of 9 quantitative lipid-related traits measured in a population-based cohort. In both examples, we compare the performance of pASTA with other competing methods as described in Table 1.

Application to Case-Control Studies of T2D

This analysis includes data from a subset of subjects in 6 independent studies of T2D: FIN-D2D 2007 (D2D2007), the DIAbetes GENetic study (DIAGEN), the Finland-United States Investigation of NIDDM Genetics Stage 2 (FUSION S2), the Nord-Trøndelag Health Study 2 (HUNT), the METabolic Syndrome In Men Study (METSIM), and the Tromsø Study (TROMSO) [25]. We investigated the associations of T2D and two SNPs in the FTO gene, rs6499640 and rs1121980, respectively, after taking into account potential interaction effects with BMI. These two SNPs were chosen as candidates because they showed significant interactions with BMI in a previous study [26]. Descriptive statistics of the 6 studies are summarized in Appendix Table A8. In particular, the total sample size of the 6 studies was 9,624, with 4,422 cases and 5,202 controls. Sample sizes of the 6 studies vary from 1,058 to 2,219, and the case to control ratios vary from 0.37 to 1.75. The minor allele frequencies of rs6499640 and rs1121980 in controls range from 0.37 to 0.43 and 0.39 to 0.48, respectively, across the studies.

The outcome of interest is the presence (D = 1) or absence (D = 0) of T2D. Each SNP is coded as G = 0, 1, 2 representing the number of minor alleles in rs6499640 or rs1121980 for each subject assuming an additive genetic susceptibility model. Let E denote BMI, A denote age in years, and S represent an indicator of sex (male = 1). For the i-th individual in the k-th study, the standard logistic regression model incorporating G-E interaction is specified as

and the standard logistic regression model of marginal SNP effect (adjusted for BMI, age, and gender) is specified as

Model 3 and Model 4 were fitted separately to the 6 studies. Then, the various meta-analysis methods were applied to the summary statistics. Note that the METSIM study only contains male subjects, so the sex indicator was removed from the two models for this study.

Application to the NFBC1966 Study

We obtained the genotype and phenotype data of the NFBC1966 study from dbGaP with accession number phs000276.v2.p1. The NFBC1966 data contain 5,402 individuals and 364,590 SNPs, with multiple metabolic traits measured on subsets of individuals. Among these traits, we considered 9 quantitative traits in the present study as phenotypes: C-reactive protein (CRP), glucose, insulin, total cholesterol (TC), high-density lipoprotein (HDL), low-density lipoprotein (LDL), triglycerides (TG), systolic blood pressure, and diastolic blood pressure. BMI was the environmental factor of interest. Following previous studies [27, 28], we excluded individuals with missing phenotype data or having discrepancies between reported sex and sex determined from the X chromosome. We excluded SNPs with a minor allele frequency <1%, having missing values in >1% of the individuals, or with a Hardy-Weinberg equilibrium p value <10^–4. This left us with 5,123 individuals and 319,147 SNPs. For each phenotype in turn, we quantile transformed the phenotypic values to a standard normal distribution, regressed out sex, oral contraceptives, and pregnancy status effects [27], and quantile transformed the residuals to a standard normal distribution again. We replaced the missing genotypes for a given SNP with its mean genotype value. To account for potential relatedness in the Finnish population, one can use a kinship matrix in a linear mixed model to account for population admixture. As an alternative, we find that using linear regression models with principal components (PCs) as covariates also effectively controls for population stratification in the data. For example, with 10 PCs, the genomic control factors for the 10 phenotypes ranged from 0.96 to 1.01, with a median value of 0.99. Therefore, we extracted the top 10 PCs from the genotype matrix as covariates to control for potential population stratifications.

The linear regression model with G-E interaction is specified as

and the linear regression model of marginal SNP effect (adjusted for BMI and PCs) is specified as

where k indexes phenotypes, i indexes subjects, E denotes BMI, G denotes the number of minor alleles, and C denotes the vector of top 10 PCs. We performed a genome-wide analysis of multiple phenotypes by fitting Model 5 and Model 6 separately to each of the 9 phenotypes for all 319,147 SNPs in the data. Then for each SNP, we calculated the meta-analytic p values using ASSET and pASTA, respectively. Correlations among the resultant Z-statistics of the 9 phenotypes were estimated using observed phenotypic correlations.

After the genome-wide screening using the 9 phenotypes separately, we performed another genome-wide association analysis by analyzing all 10 traits jointly with the multivariate linear mixed model implemented in GEMMA [29, 30] and compared the results of the two GWAS. We then selected the top 100 SNPs that showed the strongest associations for subsequent interaction analysis. We performed a SNP-BMI interaction analysis on these 100 SNPs.

Results

Simulation Study 1: Combining Signals across Multiple Independent Case-Control Studies

Type I Error

In both scenarios, empirical type I error rates at thresholds α = 10^–2, 10^–3, and 10^–5 are well-preserved for all methods except pASTA(MA+EB) (Table 2), which tends to be conservative (e.g., type I error rate being 3.8 × 10^–6 at the level of 10^–5). This has been noted in previous studies of EB-type adaptive methods [15].

Power Comparisons of ASSET and pASTA Based on 1-df Tests for Marginal Genetic Association

We extensively investigate how pASTA(MA) performs compared to ASSET when no environmental factors are considered (i.e., data generated under Model 1). The power of pASTA to detect overall genetic associations is almost the same as that of ASSET as the odds ratio of the marginal genetic effects varies, regardless of the directions of the genetic effects (Appendix Fig. A1a, A1e). However, there is an observable difference in their accuracy of determining the exact subset of non-null studies for both scenarios (Appendix Fig. A1b, A1f). When the marginal genetic effects of non-null studies exist in the same direction (in this case, positive direction), the accuracy of pASTA to determine the correct subset increases as the signal grows stronger and reaches an empirical power of 0.72 as the odds ratio goes up to 1.2. However, two-sided ASSET is rarely able to identify the exact true subset, and the corresponding power stays around 0 as the odds ratio varies. The main reason for such a low accuracy is that two-sided ASSET searches for an associated subset in both positive and negative directions, identifies a subset that gives the largest │Z(S)│ for each direction, and takes the union of the two subsets as the final identified subset regardless of the significance of association in either direction. As a result, two-sided ASSET generally identifies more false positive studies when the significant associations are all in one direction. For example, when associations only exist in the positive direction in a subset of studies, the log odds ratios in the truly null studies are likely to be estimated as negative, and thus some null studies will be selected by two-sided ASSET. When the associations exist in opposite directions, two-sided ASSET has improved power to identify the correct subset of studies, but pASTA still consistently outperforms ASSET (Appendix Fig. A1b, A1f). For both scenarios, the sensitivity of ASSET consistently stays higher than that of pASTA, while the specificity of ASSET is relatively low (Appendix Fig. A1c, A1d, A1g, A1h). These results are intuitively reasonable, since sensitivity is the proportion of non-null studies that are correctly identified and does not require that the exact subsets are identified.

When all of the 10 studies are non-null (Appendix Fig. A2a), pASTA is as powerful as ASSET in detecting signals of marginal associations, and only slightly less powerful than Fisher’s method when the odds ratio of marginal genetic effects is small (<1.1). In addition, the fixed-effect meta-analysis method outperforms the other three 1-df methods (as expected) when the marginal odds ratio is <1.1. The power becomes very similar after this effect size. When only 4 out of the 10 studies have true signals (Appendix Fig. A2b), the subset-based methods provide power gain over the fixed-effect method.

Power Comparison of All Methods with Data Generated under Model 2

Detection of Association

We consider 5 scenarios (Appendix Table A4) to evaluate the power of detecting overall genetic associations. When there exists only marginal genetic effect but no G-E interaction (Scenario 1), 1-df methods are slightly more powerful in discovering the overall associations than 2-df methods when the number of non-null studies k₀ ≤5. For example, pASTA(JOINT) reaches a power of 0.85 when 4 out of the 10 studies are non-null, which is 9.1%, 7.1%, and 9.9% lower than the power of Fisher’s(MA), ASSET, and pASTA(MA), respectively (Fig. 1a). This loss of power results from a penalty of the additional degree of freedom incorporated in the testing procedure in the absence of G-E interaction. In Scenario 2, when G-E interaction exists and the marginal genetic effect is moderate, all 2-df methods provide substantial power gain over 1-df methods (Fig. 1b). For example, when k₀ = 7, the power of pASTA(MA+EB) reaches 0.99, while that of ASSET stays around 0.6. In Scenario 3, where the marginal genetic effect is even smaller, the power of all methods is compromised, with pASTA(MA+EB) achieving the best power of 0.79 when k₀ = 7. The gold standard test has the second-best performance with a power of 0.61 when k₀ = 7. All 1-df methods barely detect any signals (Fig. 1c). When the marginal association is nearly null but subgroup effects of a genetic factor present in opposite directions in the exposed and unexposed group (Scenario 4), 2-df methods are more powerful in detecting signals than their 1-df counterparts. Specifically, pASTA(MA+EB) performs the best when k₀ > 3, already achieving a power of 0.82 when k₀ = 3. Other 2-df methods perform slightly worse, with their power ranging from 0.44 to 0.52 when k₀ = 3, but they still yield power >0.8 when k₀ > 5. All 1-df methods fail to detect associations (Fig. 1d). In a much less common situation where “cross-over” interaction exists (Scenario 5), the empirical power curves show a similar trend (Fig. 1e) to those in Scenario 4. It is instructive to note that for Scenarios 2–5, where interaction effects are present, we tested the null hypothesis of JOINT for the gold standard test. Therefore, pASTA(MA+EB) outperforms the gold standard when k₀ is large, since the EB estimator exploits the putative G-E independence and can be more efficient in this case.

Fig. 1.

Power comparison of 9 meta-analysis methods for detecting overall associations and identifying true non-null studies in case-control studies. Empirical power curves (a–e) at a significance level of 0.001 are shown for 9 meta-analysis methods (see Table 1). Power curves of identifying exact subset (f–j) are plotted only for subset-based meta-analysis method ASSET and pASTA. In each simulation, a total of 10 studies are analyzed, and each study contains 6,000 cases and 6,000 controls. Scenario 1 considers no environmental main effect or interaction, representing a situation where only marginal genetic effects exist (data generated from Model 1); Scenario 2 assumes strong interaction and moderate genetic main effect; Scenario 3 assumes strong interaction and no genetic main effect; Scenario 4 assumes positive interactions and negative genetic main effects, resulting in nearly-null marginal genetic effects; Scenario 5 considers negative interaction and positive genetic main effect.

Subset Identification

When the environmental factor is absent (Scenario 1), pASTA(MA) achieves the highest accuracy in identifying the exact subset of non-null studies, with the corresponding power being constantly around 0.4 (Fig. 1f). The sensitivities of the subset-based methods, whether they are 1-df or 2-df, are close to one another as the number of non-null studies varies, and all of the methods are able to identify >80% of the non-null studies (Appendix Fig. A3a). The specificity of ASSET is much smaller than that of the other subset-based methods based on pASTA (Appendix Fig. A3f). When G-E interaction comes into play (Scenario 2–5), pASTA(MA+EB) performs the best in identifying the exact subset and has both the highest sensitivity and specificity across the 4 scenarios. In particular, in Scenario 4 where the genetic variant has opposite effects on exposed and unexposed groups, the highest empirical probability of pASTA(MA+EB) to select the correct non-null studies is 0.46, when 3 out of the 10 studies are non-null. pASTA(JOINT) and pASTA(MA+CC) have the second-best performance in terms of determining the exact subset of non-null studies as well as sensitivity and specificity across Scenarios 2–5. It is noteworthy that in all 5 scenarios, the specificity of ASSET is the lowest among the subset-based methods, which indicates that ASSET tends to give more false positives of non-null studies than pASTA.

Simulation Study 2: A Multiple-Phenotype Analysis in a Cohort Study

Type I Error

The type I error rates are presented in Table 2. In general, the type I error rates of all methods tend to be conservative at all thresholds under evaluation. The only exception is that pASTA(JOINT) yields a type I error of 1.2×10^–5 at the level of 10^–5 when no environmental factor is considered. The conservativeness of the type I errors may be due to the inaccurate estimation of Σ, the correlation matrix of the Z-statistics.

Power Comparison

Detection of Association

Figure 2a reveals that in the absence of interaction, pASTA(MA) outperforms the other 2 methods when >5 out of the 9 traits have genetic associations. The power of ASSET drastically falls when the number of non-null traits exceeds 5 due to computational issues of the ASSET package. According to Figure 2e, when interaction exists, pASTA(JOINT) provides a great power gain in detecting the overall associations and reaches a power of 0.92 when 5 out of the 9 traits are non-null. The power of the 1-df methods stays below 0.9. We also investigate the relative power of the subset-based methods with respect to different correlation values among phenotypes: low (ρ = 0.2), moderate (ρ = 0.5), or high (ρ = 0.7). We removed ASSET from the comparison due to the non-convergence of its algorithm when phenotypic correlations are large. Appendix Figure A4a and A4d shows that pASTA has high power when the correlation between phenotypes is low (ρ = 0.2). For example, when G-E interaction exists, pASTA(JOINT) achieves a power of 0.8 when 5 out of the 9 traits are non-null (Appendix Fig. A4d). When the correlation becomes moderate (ρ = 0.5) and large (ρ = 0.7), the power of pASTA reduces, with pASTA(JOINT) achieving power around 0.60 and 0.45, respectively, in the presence of G-E interaction when 7 out of the 9 traits are non-null (Appendix Fig. A4e, A4f).

Fig. 2.

Power comparison of subset-based meta-analysis methods applied to correlated quantitative phenotypes. A cohort study with 9 correlated phenotypes measured and 6,000 subjects is simulated. The performance of the subset-based methods is evaluated in the absence (a–d) and presence (e–h) of gene-environment interaction. The level of significance is 0.001, and simulations are repeated 2,000 times.

Subset Identification

Similar to what we observe in Simulation Study 1, the accuracy of ASSET in identifying the exact subset of non-null traits is constantly close to 0 in both scenarios (Fig. 2b, f). In Figure 2c and 2g, the sensitivities of the pASTA approaches are close to one another and stay >0.8 as the number of non-null traits varies, which indicates that these methods are able to identify >80% of the non-null traits in one study. On the other hand, the sensitivity of ASSET decreases as the number of non-null traits increases. The curves of specificity (Fig. 2d, h) show similar trends to those of sensitivity.

Application 1: T2D Data

The exponential of the parameter estimates (α^(k)_G, β^(k)_G, and β^(k)_GE), their corresponding 95% confidence interval (CI), and p values obtained from Model 3 and Model 4 for each of the 6 studies are summarized in Appendix Table A9. The EB estimates of the interaction terms and their associated Wald χ² statistics were obtained using the R package CGEN [31]. For SNP rs6499640, the odds ratios of marginal genetic effect range from 0.89 to 1.11, and none of them are significant at the level of 0.05. When the interaction between SNP and BMI is included in the model, the odds ratios of genetic main effect do not change remarkably from those of marginal genetic effects. However, despite the small effect sizes of interactions in all 6 studies, the interaction effects are statistically significant in D2D2007 for the CC method (p = 0.017) and in FUSION S2 for the CC (p = 0.0006) and the EB (p = 0.003) methods. SNP rs1121980 has significant marginal associations with T2D in FUSION S2 (OR_G = 1.19, p = 0.011) and HUNT (OR_G = 1.21, p = 0.020), while none of the studies show significant SNP-BMI interactions for this SNP.

We then meta-analyzed the results from the 6 studies by applying all methods in Table 1, and we present the results of the 8 methods in Table 3. For SNP rs6499640, all 1-df methods give insignificant p values for marginal genetic associations with T2D at the level of 0.05. As for the 2-df methods, pASTA(JOINT) gives the smallest p value (p = 0.008), and pASTA(MA+CC) gives a slightly higher p value (p = 0.009). The result of pASTA(MA+EB) is not statistically significant (p = 0.090), which may be due to the presence of a G-E correlation in this example. All 2-df pASTA methods identify D2D2007 and FUSION S2 as the non-null studies, which is consistent with what we observe from the results of the 2-df tests for each study separately. We also investigated the odds ratios of T2D across different levels of BMI by fitting a logistic regression adjusting for age, gender, and the studies to which the subjects belong. Appendix Figure A5 shows that the marginal odds ratio for SNP rs6499640 is close to 1, but when G-E interaction is taken into account, this SNP reveals a significant protective effect at the level of 0.05 when BMI = 20 and 25. When BMI = 30, the odds ratio of T2D becomes positive and is significantly different from the one when BMI = 20.

Table 3.

Meta-analysis of 6 case-control studies of type 2 diabetes

When there is no evidence for G-E interactions in individual studies, as with SNP rs1121980, all 2-df methods fail to detect the overall genetic associations, but the p values obtained from 2-df pASTA are smaller than those from 2-df Fisher’s p value-combined methods. Furthermore, the p value given by pASTA(MA) (p = 0.024) is smaller than those of Fisher’s(MA) (p = 0.080) and ASSET (p = 0.031). pASTA also provides more plausible subset identification than ASSET when there are only significant marginal genetic effects. In addition to FUSION S2 and HUNT as identified by all pASTA methods, ASSET also includes D2D2007, METSIM, and TROMSO as associated studies, whose p values with respect to marginal associations in study-specific analysis (Appendix Table A9) are far from significant. This data analysis demonstrates that without sacrificing much efficiency in the absence of G-E interaction, pASTA provides a potential power gain when marginal genetic effects are modest and interactions are involved. Moreover, pASTA yields a more plausible subset of non-null studies.

Application 2: NFBC1966 Data

We performed a genome-wide multiple-phenotype analysis by applying pASTA and ASSET to the lipid-related traits in the NFBC1966 study. Appendix Table A10 presents the SNPs whose meta-analytic p values are smaller than the genome-wide significance level (p = 5 × 10^–8). Seventeen SNPs are identified by pASTA(MA) at the genome-wide significance level. ASSET identifies 7 SNPs as the genome-wide significant SNPs, and the only overlapping SNP between the two sets is rs754524, which is located in the APOB region on chromosome 2. However, it should be noted that since ASSET fails to produce p values for some of the SNPs due to computational issues, the list of genome-wide significant SNPs generated by ASSET is not comprehensive and, thus, is not comparable to the list generated by pASTA(MA). Taking G-E interaction into consideration, pASTA(JOINT) identifies 13 SNPs, which are a subset of the SNPs identified by pASTA(MA). SNP rs3764261, located on chromosome 16 near the CETP gene, is identified by pASTA(MA) (p = 3.6 × 10^–20) and pASTA(JOINT) (p = 7.1 × 10^–20) as the top SNP. Figure 3 presents the QQ plots corresponding to the 3 methods.

Fig. 3.

Quantile-quantile plots of p values obtained using ASSET, pASTA(MA), and pASTA(JOINT) in the genome-wide analysis of the North Finland Birth Cohort 1966 study. The horizontal dashed line indicates the genome-wide level of significance 5 × 10^–8.

We then compared the 17 genome-wide significant SNPs identified by pASTA(MA) to the 23 SNPs identified by another GWAS where all 10 traits were analyzed jointly using a multivariate marginal association test (i.e., GEMMA) (Appendix Table A10). Seven SNPs belong to both lists, and the top 5 SNPs identified by GEMMA and pASTA(MA), respectively, are the same (though in a slightly different order). We considered the top 100 SNPs identified by GEMMA for follow-up interaction analysis. The p values produced by GEMMA for the 100 SNPs range from 5.5 × 10^–30 to 7.0 × 10^–5. The results of the combined p values for the complete set is presented in Appendix Table A11. Two SNPs, rs2083637 (on chromosome 8, near gene LPL) and rs754524 (on chromosome 2, in the APOBregion), are used to illustrate the performance of the methods in the presence (rs2083637) or absence (rs754524) of G-E interaction. These two examples offer insight into the properties of the proposed methods and illustrate the advantage of incorporating G-E interaction into the analytic framework very clearly.

Appendix Table A12 summarizes the linear regression coefficients and their corresponding 95% CIs and p values obtained from Model 5 and Model 6. SNP rs2083637, located near the LPL gene, is marginally associated with HDL and TG, with coefficients being –0.1 (p = 3.4 × 10^–6) and 0.1 (p = 5.6 × 10^–6), respectively. CRP (p = 1.3 × 10^–3), TC (p = 6.3 × 10^–3), and LDL (p = 8.2 × 10^–3) are associated with rs2083637 through SNP-BMI interactions. On the other hand, rs754524 has a positive marginal effect on HDL (p = 3.6 × 10^–2) and a negative marginal effect on TC (p = 6.3 × 10^–10), LDL (p = 2.1 × 10^–11), and TG (p = 1.8 × 10^–2), while no significant SNP-BMI interactions are observed for these 9 phenotypes under the 0.05 threshold.

For rs2083637, ASSET and pASTA(MA) both identify TG and HDL as associated phenotypes (Table 4). However, both pASTA(JOINT) and pASTA(MA+CC) identify 3 additional phenotypes (CRP, TC, and LDL), for which rs2083637 has significant interaction with BMI, when G-E interaction is taken into account. Appendix Figure A6 presents the effect and the corresponding 95% CIs for rs2083637 on CRP, TC, and LDL when BMI equals the mean, mean ± 1 standard deviation (SD), and mean ± 2 SD. We observe variation of the SNP effect sizes across different levels of BMI. For example, rs2083637 has positive but not significant effect on CRP for subjects with average BMI, while for subjects whose BMI is 2 SD greater than the average, the association becomes negative and significant at the level of 0.05. Therefore, the incorporation of G-E interaction can enhance the knowledge that is not available in marginal association testing. The overall p value of associations given by pASTA(JOINT) (2.5 × 10^–6) is also smaller than those given by 1-df methods (2.7 × 10^–5 for pASTA[MA] and 3.6 × 10^–6 for ASSET). For rs754524, pASTA(MA) identifies 2 phenotypes (TC and LDL), both of which have significant marginal associations with this SNP in phenotype-specific analysis, as associated phenotypes. SNP rs754524 has been previously reported to be associated with both traits [32, 33]. ASSET identifies two additional phenotypes, glucose and HDL, and fails to identify TC. The evidence for the marginal association of this SNP with HDL is not strong [28], which indicates that such identified association may be false. pASTA(JOINT) identifies the same set of associated traits as pASTA(MA). This example shows that by incorporating interactions one can obtain smaller p values for detecting overall associations and discover more relevant traits that are missed by ASSET, i.e., the ones that show significant G-E interaction effects in phenotype-specific analysis.

Table 4.

Joint analysis of multiple phenotypes in the North Finland Birth Cohort 1966 study

Discussion

We propose a powerful subset-based framework for the analysis of association studies after accounting for potential G-E interactions. Simulation studies demonstrate that this framework improves the power to detect overall associations as well as the accuracy of determining associated studies/traits by incorporating G-E interactions, while maintaining comparable power to ASSET in the absence of such interactions. Our data examples exemplify that 2-df pASTA yields smaller meta-analytic p values for SNPs that show significant G-E interactions in study-/ phenotype-specific analysis compared to 1-df methods where only marginal genetic associations are considered (Tables 3 and 4). In addition, our analysis on the NFBC1966 data illustrates that pASTA is able to identify new phenotypes that are associated with given susceptibility loci after considering G-E interactions that will be missed by traditional marginal association analysis (Table 4). These properties make pASTA a powerful tool to identify candidate SNPs that are potentially associated with the trait(s) of interest in the screening step of a GWAS, which then can be followed by deeper characterization/interpretation of the associations through other methods, such as biological validation. For example, RNA sequencing experiments can be performed to examine whether the candidate SNPs affect the expression level of the target gene [34]. Crispr-Cas9 knockout screening can be carried out in human cell lines to determine the functional significance of these candidate SNPs [35].

Our simulation studies show that in the presence of interaction, pASTA(MA+EB) has the best performance in terms of all 4 evaluation metrics in a case-control setting where G-E independence holds. The rest of the 2-df methods (Fisher’s methods, pASTA[JOINT], and pASTA[MA+CC]) give similar power for detecting overall associations across the scenarios (Fig. 1b–e), but pASTA has the advantage of being able to identify the subset of studies/traits that are most likely to contribute to the associations. For a cohort study with multiple phenotypes, where G-E independence does not apply, pASTA(JOINT) gives the best performance in terms of the 4 metrics (Fig. 2e–h).

There are several limitations of our proposed approach that lend the problem to future research. The present proposal for pASTA focuses on testing the existence of genetic associations using p values as inputs and is purely bi-directional in nature. It aggregates the evidence of both genetic effects and G-E interactions from individual studies but does not specify the directions of such effects or interactions. It is a pure testing approach and does not provide pooled regression coefficient estimates from meta-analysis. Thus, it loses insight into the nature of the environmental heterogeneity across cohorts. Further work is needed to define the subsets of studies that have similar direction of genetic effects in subsets defined by the environmental factor (e.g., the sign of estimates β_G in the group E = 0 and of β_GE + β_G in the group E = 1 for a binary exposure E). The method is currently based on the assumption of fixed-effect(s) meta-analysis. Since the input is p values, the method can be extended to random-effect meta-analysis as well as hybrid methods that combine fixed- and random-effect meta-analysis [4].

In our simulation studies and data analysis, pASTA uses phenotypic correlations to approximate the correlations among Z_ks, which may result in compromised power and conservative type I error. In addition, such approximation requires individual level data on multiple traits, which are not always available in GWAS. In the case where only summary data are accessible, an alternative is to follow the idea in Liu and Lin [36] and estimate the phenotypic correlations using the summary across a large number of independent null SNPs, since the phenotypic correlations are independent of genotypes under the null for a 1-df Wald test. However, using phenotypic correlations in place of the Z-statistics correlations when our proposed Z-statistics are obtained by transforming p values of 2-df tests is only an ad hoc choice and by no means optimal. Despite the fact that our choice of using phenotypic correlations is a convenient approximation that is not theoretically proven, our simulation studies show that the properties of the resultant testing procedures are desirable. Optimizing the estimation of such correlation beyond 1-df Wald tests is an issue that needs to be further addressed.

There are some existing methods that incorporate phenotypic correlations into the construction of test statistics (as opposed to the evaluation of p values in our method). For example, Dutta et al. [37] have studied the effect of modeling kernel matrix in the score test statistic using estimated phenotype covariance for testing pleiotropic effects of rare variants based on multiple kernel regression (Multi-SKAT). While we only use phenotypic correlation for characterizing the distribution of our combined statistics instead of using it directly during the process of constructing the statistics, incorporating phenotypic correlation into the test statistics may further improve the power of our method.

We only focused on candidate loci with known marginal association or G-E interaction in our data examples. Simulation results indicate that pASTA is scalable to a genome-wide level and maintains a type I error rate at a threshold of α <10^–5. In the T2D data analysis, the average run time of pASTA for the marginal effect, i.e., pASTA(MA), for one SNP was 181.7 ms on a 2.9-GHz Intel processor, and the incorporation of G-E interaction, i.e., pASTA(JOINT), yielded a similar computation time (159.2 ms). For the data from the NFBC1966 study, where correlations across the test statistics need to be taken into account, the average run times for both pASTA(MA) and pASTA(JOINT) for one SNP were 1.7 s. For a typical GWAS to analyze one million SNPs, the projected computation time of pASTA on 20 CPU cores is approximately 1 day. ASSET was slightly faster than pASTA in the NFBC1966 data analysis (1.5 s) but twice slower in the T2D data analysis (3.3 s). Computing time was calculated using the R package microbenchmark [38]. We plan to optimize the pASTA implementation code and develop an efficient R package akin to ASSET.

Appendix

Detailed Derivation of Meta-Analytic p Value

For an observed test statistic Z_meta = T_obs, the DLM-approximated p value can be expressed as:

For a fixed subset S_γ ⊆ {1, …, K},

where h(•) is the probability distribution function of Z(S_γ). The last inequality is justified by the “separability” assumption that, given the current subset, its neighbors have independent Z-statistics.

(i) When k ∉ S_γ, given that Z(S_γ) = z,

Hence,

and

(ii) When k ∈ S_γ, given that Z(S_γ) = z,

Hence,

P(Z(S_γ – k) < z│Z(S_γ) = z) = P(Z_k >l_k(z)│Z (S_γ) = z).

This conditional probability of Z_k given (S_γ) can be evaluated based on a joint bivariate normal distribution, which is specified in the main text. For numerical evaluation of the integrals, we directly apply the R function integrate().

Metrics Used in Power Analysis

We used two main metrics for evaluation:

• Power of detecting genetic associations that truly exist, defined as the probability of rejecting the relevant null hypotheses H_0k within study k with a pre-specified level α. This quantity is estimated as the empirical proportion of replications where the meta-analytic p value combined across all studies is smaller than α;

• Power of identifying the exact true subset of non-null studies, defined as the probability of correctly identifying such a subset. The definition of non-null studies is specified above. The corresponding empirical proportion of replications declaring the correct subset to be non-null is used to estimate this power.

We also report the sensitivity and specificity of identifying the subset of non-null studies, defined as the proportion of non-null studies that are correctly identified by subset-based approaches, and the proportion of truly null studies that are declared null, respectively. They are both estimated by the corresponding empirical proportions in the simulation. Note that sensitivity and specificity allow false positives and false negatives in the identified subset, respectively, while power of identifying exactly the true subset is a much more stringent criterion.

Acknowledgement

The research of B.M. was supported by National Science Foundation grant DMS 1406712, the research of S.L. was supported by National Institute of Health grant R01-HG008773, and the research of M.B. was supported by HG009976. The FUSION study was supported by DK062370. We thank the D2D2007, DIAGEN, FUSION S2, HUNT, METSIM, and TROMSO investigators for providing access to their data. The NFBC1966 study is conducted and supported by the National Heart, Lung, and Blood Institute (NHLBI) in collaboration with the Broad Institute, UCLA, University of Oulu, and the National Institute for Health and Welfare in Finland. This manuscript was not prepared in collaboration with investigators of the NFBC1966 study and does not necessarily reflect the opinions or views of the NFBC1966 Study Investigators, Broad Institute, UCLA, University of Oulu, National Institute for Health and Welfare in Finland, and the NHLBI.

Disclosure Statement

The authors declare no conflicts of interest.

Web Resources

Current annotated code, https://github.com/youfeiyu/pASTA/.

References

---

Author Contacts

Bhramar Mukherjee

Department of Biostatistics, School of Public Health

University of Michigan, 1420 Washington Heights

Ann Arbor, MI 48109-2029 (USA)

E-Mail bhramar@umich.edu

---

Article / Publication Details

First-Page Preview

Received: May 31, 2018
Accepted: January 04, 2019
Published online: May 27, 2019
Issue release date: June 2019

Number of Print Pages: 32
Number of Figures: 3
Number of Tables: 4

ISSN: 0001-5652 (Print)
eISSN: 1423-0062 (Online)

For additional information: https://www.karger.com/HHE

---

Copyright / Drug Dosage / Disclaimer

Copyright: All rights reserved. No part of this publication may be translated into other languages, reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording, microcopying, or by any information storage and retrieval system, without permission in writing from the publisher.
Drug Dosage: The authors and the publisher have exerted every effort to ensure that drug selection and dosage set forth in this text are in accord with current recommendations and practice at the time of publication. However, in view of ongoing research, changes in government regulations, and the constant flow of information relating to drug therapy and drug reactions, the reader is urged to check the package insert for each drug for any changes in indications and dosage and for added warnings and precautions. This is particularly important when the recommended agent is a new and/or infrequently employed drug.
Disclaimer: The statements, opinions and data contained in this publication are solely those of the individual authors and contributors and not of the publishers and the editor(s). The appearance of advertisements or/and product references in the publication is not a warranty, endorsement, or approval of the products or services advertised or of their effectiveness, quality or safety. The publisher and the editor(s) disclaim responsibility for any injury to persons or property resulting from any ideas, methods, instructions or products referred to in the content or advertisements.