Gene-based association tests using GWAS summary statistics and incorporating eQTL.

Although genome-wide association studies (GWAS) have been successfully applied to a variety of complex diseases and identified many genetic variants underlying complex diseases via single marker tests, there is still a considerable heritability of complex diseases that could not be explained by GWAS. One alternative approach to overcome the missing heritability caused by genetic heterogeneity is gene-based analysis, which considers the aggregate effects of multiple genetic variants in a single test. Another alternative approach is transcriptome-wide association study (TWAS). TWAS aggregates genomic information into functionally relevant units that map to genes and their expression. TWAS is not only powerful, but can also increase the interpretability in biological mechanisms of identified trait associated genes. In this study, we propose a powerful and computationally efficient gene-based association test, called Overall. Using extended Simes procedure, Overall aggregates information from three types of traditional gene-based association tests and also incorporates expression quantitative trait locus (eQTL) information into a gene-based association test using GWAS summary statistics. We show that after a small number of replications to estimate the correlation among the integrated gene-based tests, the p values of Overall can be calculated analytically. Simulation studies show that Overall can control type I error rates very well and has higher power than the tests that we compared with. We also apply Overall to two schizophrenia GWAS summary datasets and two lipids GWAS summary datasets. The results show that this newly developed method can identify more significant genes than other methods we compared with.

Introduction

Although genome-wide association studies (GWAS) have successfully identified thousands of single nucleotide polymorphisms (SNPs) associated with a wide range of complex human traits[1,2], there is a common limitation in which GWAS focus on only a single genetic variant with a trait at a time. This limitation may limit the power to identify clinically or biologically significant genetic associations[3]. Furthermore, many genome-wide significant genetic variants are in linkage disequilibrium (LD). Different LD patterns can cause non-replicated results of the same variant in different populations[4,5]. Therefore, several powerful gene-based statistical association tests, in which the genetic information of all genetic variants in a gene is combined to obtain more informative results, have been developed, such as the Burden Test (BT)[6], the Sequence Kernel Association Test (SKAT)[7], and the Optimized SKAT (SKATO)[8].

When individual-level genotype and phenotype data are not available, the traditional gene-based association tests, BT, SKAT, and SKATO, can be extended by using GWAS summary statistics[9]. Currently, there are many GWAS summary statistics available in public resources[10]. In GWAS summary statistics, the Z-scores of genetic variants in a gene are assumed to asymptotically follow a multivariate normal distribution with a correlation matrix among all genetic variants in a gene under the null hypothesis[11], where the correlation matrix can be estimated by LD among the genetic variants in the gene[12,13]. When individual-level data are not available, LD is usually estimated using external reference panels[14,15] (i.e., 1000 Genomes Project[16]). Due to the small sample size of reference panels used to estimate LD, statistical noise (i.e., inflated type I error rates or large numbers of false positives) often exists which needs to be accounted for[17,18]. One way to reduce the statistical noise is to correct the estimated LD by a regularization procedure[19]. In the regularization procedure, a statistical white Gaussian noise is added to the LD matrix which is estimated by a reference panel. After correcting the estimated LD by the regularization procedure, we can assume that, under the null hypothesis, the Z-scores from GWAS summary statistics asymptotically follow a multivariate normal distribution with the correlation matrix being the corrected LD matrix among the genetic variants in a gene.

To increase statistical power in identifying genes that are associated with complex diseases, PrediXcan[20] and transcriptome-wide association study[12,21] (TWAS) were developed by incorporating expression quantitative trait locus (eQTL) data into GWAS. As pointed out by Zhang et al.[15], PrediXcan and TWAS can be viewed as a simple weighted linear combination of genetic variants with an eQTL—derived weight. In fact, the genetic architecture of complex traits is rarely known in advance and is likely to vary from one region to another across the genome and from one trait to another[15]. Therefore, only considering one single eQTL—derived weight, such as in PrediXcan and TWAS, may lose statistical power in identifying significant genes. Zhang et al.[15] developed an omnibus test (OT) using Cauchy combination method to integrate association evidence obtained by BT, SKAT, and SKATO based on GWAS summary data with multiple eQTL‐derived weights. They showed that OT using multiple eQTL—derived weights had some potential advantages.

Inspired by the advantage of OT, in this paper, we propose a more powerful and computationally efficient method, called Overall, to aggregate the information from three types of traditional gene-based association tests (BT, SKAT, SKATO) with multiple eQTL—derived weights using GWAS summary statistics. To combine information from the three gene-based association tests, the Overall method utilizes the extended Simes procedure[5,22]. To apply the Overall method, we first need to estimate the correlation matrix among the three gene-based association tests with eQTL—derived weights under the null hypothesis. We provide an estimation method using a replication procedure[23,24]. The replication procedure only needs to be performed once to obtain the correlation matrix for each gene. Then, the p-values of Overall can be analytically computed without using permutations. To calculate the p-values of the three types of gene-based association tests (BT, SKAT, SKATO) using GWAS summary statistics with eQTL—derived weights, we use the “sumFREGAT” package in R (https://cran.r-project.org/web/packages/sumFREGAT/index.html)[9]. Once we obtain the p-values of these three tests, the p-value of our proposed method can be easily calculated using its theoretical distribution. Extensive simulation studies show that Overall can control type I error rates well and has higher power than the comparison methods across most of the simulation settings. Similar to Zhang et al.[15], we apply our method to two schizophrenia (SCZ) and two lipids trait (HDL) GWAS summary data sets. Compared with OT and other tests, the proposed method can identify more significant genes. More interestingly, some significant genes reported by GWAS catalog are only identified by our proposed method.

Statistical models and methods

Statistical models

Consider a set of genetic variants in a gene. Let be an vector of Z-scores of the genetic variants. Note that the Z-scores is either directly provided by publicly available GWAS summary statistics or calculated from a GWAS individual-level genotype and phenotype data set. We are interested in testing the null hypothesis that none of the genetic variants in the gene is associated with a trait, whereas the alternative hypothesis is that at least one genetic variant in the gene is associated with a trait. Following Svishcheva et al.[9], Gusev et al.[12], and Yang et al.[25], we assume under the null hypothesis, where is the correlation matrix among , which can be estimated by LD among the genetic variants in the gene[12,13]. If individual-level data are not available, LD can be estimated using external reference panels (i.e., 1000 Genomes Project[16]). However, if the sample size of a reference panel is small, LD may not be estimated correctly so that it will induce statistical noise (i.e., inflated type I error rates or large numbers of false positives)[17,18]. One way to correct the estimated LD is to use a regularization procedure by adding a statistical white Gaussian noise[9,19]. Let be an identity matrix, and the corrected correlation matrix can be defined aswhere is a scalar tuning parameter which represents the coefficient of proportionality between the corrected correlation matrix and the original estimated using an external reference panel. The optimal tuning parameter can be estimated by maximizing the log-likelihood function of the distribution of , that is,

Then the corrected correlation matrix . Therefore, under the null hypothesis, we consider .

Suppose that there are a total of different eQTL—derived weights from gene expression data (i.e., Genotype-Tissue Expression (GTEx) project (https://gtexportal.org/home/)), denoted as for , where represents a status without using any weights. In order to avoid the influence of the scale among genetic variants within each weight, we first standardize the eQTL—derived weights as for . Based on the th standardized weight , the weighted Z-score follows a multivariate normal distribution. That is,

We extend the three types of gene-based association tests, BT[6], SKAT[7], and SKATO[8], to incorporate the eQTL—derived weights based on GWAS summary statistics[9,26]. For the kth eQTL—derived weight, the three gene-based test statistics can be written aswhere is an vector with elements of all 1s. Under the null hypothesis, follows a distribution with 1 degree of freedom; follows a weighted sum of distributions with 1 degree of freedom; and follows a mixture of distribution[8]. The p-values of these three test statistics can be easily calculated using the “sumFREGAT” package in R (https://cran.r-project.org/web/packages/sumFREGAT/index.html)[9].

Overall method

To aggregate information from these three gene-based association tests with multiple eQTL—derived weights, we develop a novel method, called Overall, which utilizes the extended Simes procedure[5,22]. Let be the p-values of BT, SKAT, SKATO with kth eQTL—derived weight, , respectively, where denotes a status without using any weights. Thus, there are a total of p-values from these three gene-based association tests with different weights. Let be a sequence of the ascending p-values, where and . Overall combines these p-values using the extended Simes procedure[5,22], and the p-value of Overall is defined aswhere is the effective number of p-values among the gene-based association tests with multiple weights, is the th element of the ascending p-values, and is the effective number of p-values among the top association tests. We use a more robust measure to obtain the effective numbers and , which was proposed by Li et al.[5]. The values of and can be estimated as where denotes the ith eigenvalue of the correlation matrix of p-values from association tests with multiple weights (the estimation of will be discussed in the next section), is an indicator function. If the association tests are independent, all eigenvalues equal 1, and for ; if the association tests are perfectly dependent, then which is the number of tests used to calculate and the other eigenvalues equal 0. In this case, for .

The R codes and a sample data set for the implementation of Overall are available at github https://github.com/xueweic/Overall.

Estimation of under the null hypothesis

To apply our proposed method, we need to estimate the correlation matrix of p-values under the null hypothesis. Since the exact correlations among all gene-based association tests are unknown, we perform the estimation procedure with replications. For each replicate , , we implement the following two steps:

Step 1: We first generate a new Z-score vector under the null hypothesis. That is, follows a multivariate normal distribution with mean and variance–covariance matrix , where can be estimated by LD among the genetic variants in a gene using external reference panels (i.e., 1000 Genomes Project).

Step 2: We use the regularization procedure to obtain the corrected correlation matrix of Z-scores . Then, we calculate and the corresponding p-values using the simulated for . The distributions of depend on the corrected correlation matrix , and the standardized eQTL—derived weights for .

To estimate the correlation matrix of p-values used in the Overall method, we use the sample correlation matrix of the p-values obtained from the replications. We denote the sample correlation matrix of p-values as . For example, is the (1,2)-element of which is the estimated correlation between BT and SKAT without using any weights. If we let be a vector of the p-values of BT without using any weights and be a vector of the p-values of SKAT without using any weights obtained from the replications, then the sample correlation of p-values between these two tests is defined as , where is the sample correlation.

The estimation procedure to estimate is independent of our proposed method, therefore we only need to perform this procedure once for each gene. After we estimate , the p-value of Overall can be computed analytically without using permutations.

Simulation studies

Materials and comparison methods

In our studies, we use four data sets to obtain the eQTL—derived weights downloaded from the functional summary-based imputation website (http://gusevlab.org/projects/fusion/#reference-functionaldata). The resources to obtain the four eQTL—derived weights are listed in Table 1. For each eQTL data set, we use the weights estimated by the Best Linear Unbiased Prediction (BLUP)[27].

Table 1

Resources of the four eQTL—derived weights used in the simulation studies.

Study	Tissue	# of samples	References
NTR	Peripheral blood	1247	Wright et al.[28]
YFS	Whole blood	1264	Gusev et al.[12]
METSIM	Adipose	563	Gusev et al.[12]
CMC	Brain	452	Gusev et al.[12]

We compare our proposed method with three existing methods, OT[15], S-PrediXcan[29], and S-TWAS[12]. These three methods are all based on GWAS summary statistics and incorporate eQTL‐derived weights. Here, we briefly introduce these three methods.

OT: For a total of different eQTL—derived weights and the three gene-based association tests (BT, SKAT, SKATO), OT aggregates information from different weights and tests by using the Cauchy combination method[30]. The test statistic of OT is defined as and the corresponding p-value of the test statistic can be approximated by .

S-PrediXcan: For a given eQTL‐derived weight, provided by a matrix , the test statistic of S-PrediXcan is defined as , where is the estimated standard deviation of the SNP in a gene and is the estimated standard deviation of the predicted expression of a gene. The p-value of S-PrediXcan can be computed as , where is the standard normal CDF function.

S-TWAS: For a given eQTL‐derived weight, provided by a vector , the test statistic of S-TWAS is defined as , where is the estimated LD structure among the genetic variants in a gene and the corresponding p-value can be calculated by .

The number of replications needed in estimation of

To apply our proposed method, we first need to estimate the correlation matrix of p-values, , under the null hypothesis for each gene. Following the estimation procedure introduced in the method section, we generate Z-scores instead of generating individual-level genotype and phenotype data. To determine the number of replications needed in the estimation of , we consider 18 genes that contain different numbers of SNPs and have different LD structures. Supplementary Table S1 gives a summary of these 18 genes. We can see from Supplementary Table S1, the number of SNPs in a gene is ranging from 23 to 359 and the average per-SNP LD score in a gene is ranging from 12.72 to 170.85. We simulate a Z-score vector from a multivariate normal distribution with mean and variance–covariance matrix , , where is the LD matrix of each gene which can be estimated using the 1000 Genomes Project (unrelated Europeans in 1000 Genomes in Phase 3; ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/). First, we use replications to estimate under the null hypothesis, where the estimated matrix is denoted by . Then, we denote as the correlation matrix of p-values by using replications. We vary the value of from 16 to 5000, and test the null hypothesis that the two correlation matrices, and , are the same by using “lavaan” package (https://CRAN.R-project.org/package=lavaan)[31]. Supplementary Figure S1 shows that the p-values for the hypothesis testing in each gene are greater than 0.05 after replications for all of the 18 genes. Therefore, we recommend using 1000 replications to obtain for each gene under the null hypothesis. Consequently, 1000 replications are used in the following sessions to evaluate the type I error rates and powers of Overall.

Type I error rates

To evaluate if our proposed method can control type I error rates, we perform simulations based on the aforementioned 18 genes. For each of the 18 genes, we generate Z-score vectors under the null hypothesis, , where is the LD matrix of the gene estimated using the 1000 Genomes project. Then, we use the regularization procedure to obtain the corrected correlation matrix of Z-scores , and calculate the three types of gene-based association tests, BT, SKAT, and SKATO, with or without the four eQTL—derived weights (NTR, YFS, METSIM, CMC) based on the corrected correlation matrix . Finally, we apply our proposed method to combine the p-values using the estimated correlation matrix of p-values, , with 1000 replications.

We generate simulated data to mimic real lipids data which we will use in “Real data analysis” section. Gene AGTRAP is associated with lipids trait HDL[15], There are a total of 23 genetic variants in gene AGTRAP. The LD block structure of these 23 genetic variants is shown in Supplementary Fig. S2. Supplementary Figure S3 shows the estimated correlation matrix for this gene. We use replications to evaluate type I error rates of Overall for gene AGTRAP at , , , ,, and significance levels. With replications, a Bonferroni corrected significance level of can be reached to obtain the empirical type I error rates (i.e., for 28,625 genes in the real data analysis section, the Bonferroni corrected significance level is at 5% significance level). We further evaluate type I rates based on the other 17 genes. To save computational time, we use replications to evaluate type I error rates of Overall for the 17 genes at significance levels of , , and . Table 2 and Supplementary Table S2 show the estimated type I error rates of Overall under various nominal significance levels for gene AGTRAP and the other 17 genes, respectively. From these tables, we can see that our proposed method can control type I error rates very well at different significant levels.

Table 2

Estimated type I error rates at different significance levels with replications. The subscript denotes BT, SKAT, and SKATO using eQTL—derived weights; CMC, METSIM, NTR, and YFS indicate the resources to obtain the eQTL—derived weights. 0 indicates the methods without using eQTL—derived weights.

α-level	5×10^-2	1×10^-2	1×10^-3	1×10^-4	1×10^-5	1.75×10^-6
BT0	5.03×10^-2	1.06×10^-2	1.00×10^-3	1.01×10^-4	9.76×10^-6	1.84×10^-6
SKAT0	5.24×10^-2	1.07×10^-2	1.01×10^-3	1.00×10^-4	1.04×10^-5	1.80×10^-6
SKATO0	4.58×10^-2	9.57×10^-3	1.02×10^-3	1.04×10^-4	9.72×10^-6	1.46×10^-6
BTCMC	5.17×10^-2	1.04×10^-2	1.01×10^-3	9.82×10^-5	9.58×10^-6	1.72×10^-6
SKATCMC	5.08×10^-2	9.89×10^-3	9.71×10^-4	9.75×10^-5	9.48×10^-6	1.66×10^-6
SKATOCMC	5.16×10^-2	1.09×10^-2	1.17×10^-3	1.21×10^-4	1.22×10^-5	2.14×10^-6
BTMETSIM	5.02×10^-2	1.03×10^-2	1.02×10^-3	1.01×10^-4	9.86×10^-6	1.66×10^-6
SKATMETSIM	5.30×10^-2	1.08×10^-2	1.02×10^-3	9.91×10^-5	1.00×10^-5	2.12×10^-6
SKATOMETSIM	4.84×10^-2	1.05×10^-2	1.11×10^-3	1.09×10^-4	1.06×10^-5	1.84×10^-6
BTNTR	5.02×10^-2	1.06×10^-2	1.00×10^-3	9.93×10^-5	1.01×10^-5	1.76×10^-6
SKATNTR	5.09×10^-2	1.03×10^-2	9.98×10^-4	1.00×10^-4	1.01×10^-5	2.00×10^-6
SKATONTR	5.08×10^-2	1.18×10^-2	1.34×10^-3	1.45×10^-4	1.52×10^-5	2.92×10^-6
BTYFS	5.10×10^-2	1.02×10^-2	9.95×10^-4	9.95×10^-5	1.05×10^-5	2.10×10^-6
SKATYFS	4.98×10^-2	1.03×10^-2	9.97×10^-4	1.01×10^-4	1.02×10^-5	2.06×10^-6
SKATOYFS	5.58×10^-2	1.32×10^-2	1.43×10^-3	1.55×10^-4	1.69×10^-5	3.50×10^-6
Overall	4.67×10^-2	1.01×10^-2	1.12×10^-3	1.14×10^-4	1.24×10^-5	2.44×10^-6

Power comparison

To evaluate the performance of the Overall method, we use several simulations to compare the power of Overall with the power of OT, S-PrediXcan, S-TWAS, and three types of gene-based association tests with and without eQTL—derived weights. We use BEST to represent the test with the maximum power among the three traditional gene-based association tests with and without an eQTL—derived weight, S-TWAS.B and S-PrediXcan.B to represent the maximum power of S-TWAS and S-PrediXcan with each of the eQTL—derived weights, respectively. Following the simulation settings in Nagpal et al.[32] and Zhang et al.[15], we generate individual-level genotypes, phenotypes, and different gene expression levels using the following steps:

The genotype data are generated using the haplotypes of a gene obtained from the 1000 Genomes Project reference panel. To generate the genotype of an individual, , we select two haplotypes according to the haplotype frequencies from the haplotype pool and then remove genetic variants with MAF < 0.05.

We consider different weights derived from gene expression data which can be estimated using BLUP. To generate a vector of weights, , for the kth gene expression level, we randomly select causal variants according to the proportion of causal variants, . Then, the effect sizes for the kth gene expression levels and causal variants can be generated from a standard normal distribution, for , where ; otherwise, . After we rescaled the weights to ensure the targeted expression heritability , we generate the kth gene expression level by with each element of random error follows .

Let be the matrix of gene expression levels. Phenotypes are generated by using a formula with each element of random error follows , where is a vector of genetic effect sizes which can be assigned based on the phenotypic heritability .

The Z-score vector is estimated from individual-level genotype and phenotype data using beta coefficient and its standard deviation estimated based on the ordinary least squares method in linear regression.

In our simulation studies for power comparison, we consider two genes, AGTRAP and C3orf22, from the 18 genes used in the type I error evaluation and and eQTL—derived weights. AGTRAP contains 458 haplotypes for 23 genetic variants (11 common variants and 12 rare variants; MAF ranging from 0 to 0.39775); C3orf22 contains 295 haplotypes for 42 variants (18 common variants and 24 rare variants; MAF ranging from 0 to 0.43558). Supplementary Figure S2 shows the LD block structure of the 23 genetic variants at AGTRAP and the 42 genetic variants at C3orf22. We vary the proportion of causal variants with values for AGTRAP and for C3orf22. We also consider two different directions of genetic effects: (Scenario 1: Uni-directional effects) and (Scenario 2: Bi-directional effects). For each simulation scenario, we vary the proportion of gene expression heritability and the phenotypic heritability with different values of and . We consider the sample size to be 2000 (unless it is specified) and the power is calculated as the proportion of 1000 replications with p-value .

Figure 1 (Supplementary Fig. S4) show the power comparisons based on gene AGTRAP (and C3orf22) with under the Uni-directional effects () with different . We consider two settings here. First, we vary phenotypic heritability with a fixed expression heritability (Fig. 1a and Supplementary Fig. S4a). Second, we vary the expression heritability with a fixed phenotypic heritability (Fig. 1b and Supplementary Fig. S4b). Figure 2 (Supplementary Fig. S5) shows power comparisons based on gene AGTRAP (and C3orf22) under the Bi-directional effects () with different for . We also consider two simulation settings, power against the phenotypic heritability with a fixed expression heritability and power against the expression heritability with a fixed phenotypic heritability . The pattern of the power in Fig. 2 (Supplementary Fig. S5) is similar to what we observe in Fig. 1 (Supplementary Fig. S4). These figures show that (1) Overall and OT perform uniformly better than BEST, S-TWAS.B, and S-PrediXcan.B. We can see that Overall and OT boost power significantly due to integrating association evidence by different traditional tests and multiple eQTL—derived weights. Overall is slightly more powerful than OT in all of the scenarios. (2) Among BEST, S-TWAS.B, and S-PrediXcan.B, BEST is more powerful than S-TWAS.B and S-PrediXcan.B in all of the scenarios for gene C3orf22; For gene AGTRAP, S-TWAS.B and S-PrediXcan.B perform better than BEST when the proportion of causal variants in a gene is small (); otherwise, BEST performs better than S-TWAS.B and S-PrediXcan.B.

Figure 1

Power comparisons of gene-based association tests at significance level under Uni-directional effects () with based on gene AGTRAP. (a) Estimated power against phenotypic heritability with fixed expression heritability ; (b) Estimated power against expression heritability with fixed phenotypic heritability .

Figure 2

Power comparisons of gene-based association tests at significance level under Bi-directional effects () with based on gene AGTRAP. (a) Estimated power against phenotypic heritability with expression heritability ; (b) Estimated power against expression heritability with phenotypic heritability .

To evaluate if Overall and OT that integrate different types of association tests and multiple eQTL—derived weights are robust for more eQTL studies, we also consider 20 () eQTL—derived weights under Uni-directional effect and Bi-directional effect models on gene C3orf22 with settings similar to the settings in Supplementary Figs. S4 and S5. After integrating traditional gene-based association tests, we observe that the patterns of the power for are similar to that in Supplementary Figs. S4 and S5 with , and the power gain of Overall and OT is higher than that of the tests only consider one eQTL—derived weight, such as BEST, S-PrediXcan.B, and S-TWAS.B (Supplementary Fig. S6).

Furthermore, we consider simulation settings with noise to the eQTL. We consider simulation settings by adding less noise to the eQTL from the most relevant tissues and more noise to those from the less relevant tissues. For the Uni-direction scenario, we consider the first study being the most relevant tissue, where and ; depends on the phenotypic heritability . For the Bi-direction scenario, we consider the first and third studies being the most relevant tissues that have opposite effect directions, where and . Other parameter settings are the same as these in Supplementary Figs. S4 and S5. The power comparison results are shown in Supplementary Figs. S7 and S8. From these figures, we find that the patterns of the power in Supplementary Figs. S7 and S8 are very similar to those in Supplementary Figs. S4 and S5.

In all of the previous power comparisons, we use a sample size of 2000. We also consider simulation settings as those in Supplementary Figs. S7 and S8, but with a large sample size of 100,000. Supplementary Figure S9 shows the results of power comparisons. We can see from this figure, all powers are increased with this larger sample size, but the patterns of the power are very similar to those in Supplementary Figs. S7 and S8.

To remove noise in LD matrix computed from a reference sample, we shrink the observed LD matrix toward an identity matrix with the shrinkage parameter estimated by maximum likelihood. To evaluate how well this regulation process performs, we compare the powers of three traditional gene-based association tests with and without eQTL—derived weights, OT, and Overall based on corrected and uncorrected LD structure. We use the same simulation settings as those in Supplementary Figs. S7 and S8. Supplementary Figure S10 shows the power comparison results based on gene C3orf22 under Uni-directional effects and Bi-directional effects with noise to eQTL. We can see that the powers of these tests based on corrected LD structure perform better than those based on uncorrected LD structure in most of the settings.

Real data analysis

To evaluate the performance of our proposed method, we apply Overall, OT, the three traditional tests with and without eQTL—derived weights, S-PrediXcan, and S-TWAS to the GWAS summary statistics data sets used in Zhang et al.[15]: two SCZ GWAS summary data sets and two lipid GWAS summary data sets. We estimate the LD between genetic variants using the 1000 Genomes Project reference panel[16], and obtain the corrected matrix of Z-score after the regularization procedure. We consider four eQTL—derived weights estimated by the BLUP method using the resources listed in Table 1 (NTR, YFS, METSIM, CMC).

Application to the SCZ GWAS summary data

We consider two SCZ GWAS summary data sets, SCZ1 and SCZ2, which can be downloaded from the Psychiatric Genomics Consortium website (https://www.med.unc.edu/pgc/results‐and‐downloads/)[33]. SCZ1 is a meta-analysis of SCZ GWAS data set with 13,833 cases and 18,310 controls. SCZ2 is a more recent and larger SCZ GWAS summary data set with 36,989 cases and 113,075 controls for partial validation[34]. In our real data analysis, we define a gene to include all of the SNPs from 20 kb upstream to 20 kb downstream of the gene and test the association between each gene and the trait. We consider all genes according to the GENCODE version 35 (GRCh37) human comprehensive gene annotation list which can be downloaded from the GENCODE website (https://www.gencodegenes.org/human/release_35lift37.html).

To make fair comparisons among all these weighted tests, the genetic variants are removed if there is at least one weight missing in the four eQTL—derived weights. After pruning, there are 26,575 genes in SCZ1 and 17,823 genes in SCZ2 left in our final analyses. Therefore, the Bonferroni corrected significance level for gene-based association analysis is defined as 0.05 divided by the number of genes. First, we apply BT, SKAT, and SKATO with and without an eQTL—derived weight, OT, Overall, S-PrediXcan, and S-TWAS to the SCZ1 and SCZ2 data sets. Table 3 (SCZ1 and SCZ2) shows the number of genes identified by each method for the SCZ data sets, respectively. As we can see in Table 3, Overall identifies more genes than all of the other methods for two SCZ GWAS summary data sets. Among the three types of gene-based association tests, BT, SKAT, and SKATO, with or without different eQTL—derived weights, SKATO0 identifies the largest number of genes. S-TWASYFS and S-PrediXcanYFS identify the largest number of genes compared with S-TWAS and PrediXcan based on the other three eQTL—derived weights, respectively. Therefore, in Fig. 3, we only show the number of genes identified by Overall, OT, SKATO0, S-PrediXcanYFS, and S-TWASYFS. The number below each method indicates the total number of genes identified by the corresponding method. From Fig. 3, we can see that Overall identifies all of the genes identified by OT for SCZ1; for SCZ2, there are two genes identified by OT but failed to be identified by Overall; there are 66 and 24 genes identified only by Overall for SCZ1 data and SCZ2, respectively.

Table 3

The numbers of genes identified by each method for the two SCZ data sets. The subscript denotes BT, SKAT, and SKATO using eQTL—derived weights; CMC, METSIM, NTR, and YFS indicate the resources to obtain the eQTL—derived weights. 0 indicates the methods without using any weights.

	SCZ1	SCZ2	SCZoverlap	GWASSCZ1	GWASSCZ2
BT0	97	166	7	1	38
SKAT0	47	305	20	15	153
SKATO0	136	394	27	15	153
BTCMC	44	137	2	1	56
SKATCMC	12	225	6	1	134
SKATOCMC	30	263	2	1	130
BTMETSIM	44	136	5	1	48
SKATMETSIM	23	223	9	4	132
SKATOMETSIM	31	205	3	0	100
BTNTR	48	119	7	6	48
SKATNTR	27	230	9	8	141
SKATONTR	40	280	8	6	143
BTYFS	89	166	14	1	53
SKATYFS	20	223	6	7	137
SKATOYFS	47	321	7	0	140
S-PrediXcanCMC	42	43	7	0	38
S-PrediXcanMETSIM	41	44	8	1	30
S-PrediXcanNTR	48	70	14	6	59
S-PrediXcanYFS	83	128	29	2	72
S-TWASCMC	33	45	6	0	43
S-TWASMETSIM	36	29	5	1	20
S-TWASNTR	37	54	13	6	46
S-TWASYFS	64	105	29	2	58
OT	133	522	17	6	166
Overall	271	559	45	16	167

SCZ1 indicates the number of genes identified by each method for SCZ1 data; SCZ2 indicates the number of genes identified by each method for SCZ2 data; SCZoverall indicates the number of overlapping genes identified by both SCZ1 and SCZ2 data sets; GWASSCZ1 and GWASSCZ2 indicate the numbers of genome-wide significant genes that are reported in the GWAS catalog and are also identified by each method for SCZ1 and SCZ2, respectively.

Figure 3

Venn diagram of the number of genes identified by Overall, OT, and SKATO0, S-PrediXcanYFS, and S-TWASYFS for SCZ1 data (left) and SCZ2 data (right). The number below each method indicates the total number of significant genes identified by the corresponding method.

We further investigate the 90 genes identified only by Overall for the SCZ data sets by searching the GWAS catalog (https://www.ebi.ac.uk/gwas/). Among the 66 genes for the SCZ1 data set, there are six genes reported in the GWAS catalog; among the 24 genes for the SCZ2 data set, there are six genes reported in the GWAS catalog (Table 4). We also use these two SCZ GWAS data sets for partial validation. Table 3 shows that there are 45 overlapping genes identified by Overall using SCZ1 and SCZ2 data sets and only 17 overlapping genes identified by OT using both SCZ1 and SCZ2 data sets. Furthermore, we search for genome-wide significant SNPs () from the two SCZ GWAS summary data sets and consider the genes covering at least one genome-wide significant SNP from 20 kb upstream to 20 kb downstream of the gene. There are 63 genome-wide significant genes for SCZ1, and 2422 genome-wide significant genes in SCZ2. Table 3 (GWASSCZ1 and GWASSCZ2) summarizes the numbers of genome-wide significant genes that are identified by each method for the two SCZ data sets. Among the 63 genome-wide significant genes for the SCZ1 data set, Overall identifies the largest number of genes, followed by SKAT0 and SKATO0; OT, S-PrediXcanNTR and S-TWASNTR only identify 6 genes. Meanwhile, among 2422 genome-wide significant genes for SCZ2, Overall identifies 167 genes; OT identifies 166 genes; SKATO and SKATO0 identify 153 genes; S-TWASYFS and S-PrediXcanYFS only identify 58 and 72 genes respectively.

Table 4

Genes identified only by Overall based on the two SCZ data sets that are reported in the GWAS catalog.

Gene	Data	Overall	References
RAI1	SCZ1	2.63E−31	Pardiñas et al.[35]
SLC7A6	SCZ1	2.17E−15	Ikeda et al.[36]; Li et al.[37]
AP001931.2	SCZ1	1.27E−13	Schizophrenia Working Group of the Psychiatric Genomics Consortium[34]; Goes et al.[38]; Ikeda et al.[36]; Li et al.[37]; Lam et al.[39]; Periyasamy et al.[40]; Lee et al.[41]; The Autism Spectrum Disorders Working Group of the Psychiatric Genomics Consortium[42]; Pardiñas et al.[35]
MARK2	SCZ1	2.64E−07	Goes et al.[38]
GULOP	SCZ1	1.24E−07	Pardiñas et al.[35]; Ikeda et al.[36]; Li et al.[37]; Goes et al.[38]; Lam et al.[43]
ZBED4	SCZ1	9.02E−07	Goes et al.[38]
RAB11FIP5	SCZ2	1.05E−06	Goes et al.[38]; Lam et al.[43]
AL669918.1	SCZ2	2.03E−06	Goes et al.[38]
YPEL1	SCZ2	2.80E−06	Goes et al.[38]
LINC00606	SCZ2	2.57E−06	Goes et al.[38]
ERLIN1	SCZ2	2.34E−06	Goes et al.[38]
AC024597.1	SCZ2	2.56E−06	Lam et al.[39]

Application to the lipids GWAS summary data

We consider two lipids GWAS summary data sets, HDL1 and HDL2, which can be downloaded at the Center for Statistical Genetics (CSG) at the University of Michigan. HDL1 is a meta-analysis of HDL GWAS data set with about 100,000 samples downloaded at the website (http://csg.sph.umich.edu/willer/public/lipids2010/)[44]. HDL2 is the follow-up data with about 189,000 samples for partial validation downloaded at the Global Lipids Genetics Consortium (http://csg.sph.umich.edu/willer/public/lipids2013/)[45]. We perform the same analysis as we did in the previous section for the two SCZ GWAS summary data sets. After pruning and removing the genetic variants with missing weights, there are 17,389 genes in HDL1 and 16,917 genes in HDL2. Table 5 (HDL1 and HDL2) shows the number of genes identified by each method for the two lipids data sets, respectively. As we can see from Table 5, among the three traditional gene-based association tests with and without eQTL—derived weights, SKATO0 and BT0 identify the largest number of genes in HDL1 and HDL2, respectively; Among the four S-PrediXcan tests, S-PrediXcanYFS and S-PrediXcanCMC identify the largest number of genes in HDL1 and HDL2, respectively; for the four S-TWAS tests, S-TWASYFS and S-TWASCMC identify the largest number of genes in HDL1 and HDL2, respectively. For the HDL1 data set, Overall identifies the largest number of genes (249), followed by OT that identifies 233 genes; for the HDL2 data set, BT0 identifies the largest number of genes (836), followed by Overall and OT, where Overall identifies 765 genes and OT identifies 688 genes. In Fig. 4, we compare genes identified by SKATO0, S-PrediXcanYFS, and S-TWASYFS, along with Overall and OT for the HDL1 data set and genes identified by BT0, S-PrediXcanCMC, S-TWASCMC, Overall, and OT for the HDL2 data set. Again, we observe that Overall identifies the largest number of genes for the HDL1 data set and the second most for the HDL2 data set; all genes identified by OT are also identified by Overall; 82 and 24 genes are identified only by Overall and OT for the HDL1 and HDL2 data sets, respectively; there are 13 and 6 genes only identified by Overall for the HDL1 and HDL2 data sets, respectively. We search the GWAS catalog (https://www.ebi.ac.uk/gwas/). Table 6 shows that five out of 13 genes identified only by Overall based on HDL1 data have been reported, and one out of 6 genes has been reported on HDL2 data in the GWAS catalog. We also use these two HDL GWAS data sets for partial validation by looking for the number of overlapping genes identified by both of the data sets (Table 5, HDLoverlap). There are 177 overlapping genes identified by Overall for both SCZ1 and SCZ2 data sets and 167 overlapping genes identified by OT for both SCZ1 and SCZ2 data sets.

Table 5

The number of genes identified by each method for the two lipids data sets. The subscript denotes BT, SKAT, and SKATO using eQTL—derived weights; CMC, METSIM, NTR, and YFS indicate the resources to obtain the eQTL—derived weights. 0 indicates the methods without using any weights.

	HDL1	HDL2	HDLoverlap	GWASHDL1	GWASHDL2
BT0	95	836	78	50	185
SKAT0	116	174	114	99	157
SKATO0	157	762	138	104	190
BTCMC	79	130	41	46	107
SKATCMC	105	159	99	95	146
SKATOCMC	130	177	103	96	150
BTMETSIM	83	160	59	58	111
SKATMETSIM	120	259	118	102	149
SKATOMETSIM	131	199	118	98	152
BTNTR	78	136	50	49	111
SKATNTR	105	156	100	90	148
SKATONTR	131	183	111	95	154
BTYFS	88	154	50	53	113
SKATYFS	106	148	102	94	137
SKATOYFS	142	185	112	99	144
S-PrediXcanCMC	43	213	18	29	114
S-PrediXcanMETSIM	45	201	23	30	118
S-PrediXcanNTR	33	187	14	19	108
S-PrediXcanYFS	69	195	25	31	117
S-TWASCMC	40	207	17	23	109
S-TWASMETSIM	37	202	16	15	112
S-TWASNTR	25	176	10	11	97
S-TWASYFS	59	183	24	29	115
OT	233	688	167	120	190
Overall	249	765	177	122	192

HDL1 indicates the number of genes identified by each method for HDL1 data; HDL2 indicates the number of genes identified by each method for HDL2 data; HDLoverall indicates the number of overlapping genes identified by both HDL1 and HDL2 data sets; GWASHDL1 and GWASHDL2 indicate the numbers of genome-wide significant genes that are reported in the GWAS catalog and are also identified by each method for HDL1 and HDL2, respectively.

Figure 4

Venn diagram of the number of genes identified by Overall, OT, SKATO0, S-PrediXcanYFS, and S-TWASYFS for HDL1 data (left) and Overall, OT, BT0, S-PrediXcanCMC, and S-TWASCMC for HDL2 data (right). The number below each method indicates the total number of significant genes identified by the corresponding method.

Table 6

Genes identified only by Overall based on the two lipids data sets that are reported in the GWAS catalog.

Gene	Data	Overall	References
AP002954.1	HDL1	2.27E−11	Emilsson et al.[46]
EDC4	HDL1	1.65E−11	Lettre et al.[47], Kilpeläinen et al.[48], Wojcik et al.[49]
PACSIN1	HDL1	2.24E−06	Liu et al.[50]
AFF1	HDL1	2.10E−06	Spracklen et al.[51], De Vries et al.[52], Hoffmann et al.[53], Ripatti et al.[54], Richardson et al.[55]
AC106779.1	HDL1	2.85E−06	Noordam et al.[56]
NHLRC2	HDL2	1.98E−06	Hoffmann et al.[53], Richardson et al.[55], Klarin et al.[57], Qi et al.[58], Klimentidis et al.[59], Liu et al.[60]

Same as the analyses for the SCZ GWAS summary data sets, we search for genome-wide significant SNPs () from the two lipids GWAS summary statistics. There are 1911 genome-wide significant genes for HDL1 and 2682 genome-wide significant genes for HDL2. Table 5 (GWASHDL1 and GWASHDL2) summarizes the numbers of genome-wide significant genes that are identified by each method for the two lipids data sets. Among the 1911 genome-wide significant genes for the HDL1 data set, Overall identifies the largest number of genes (122), followed by OT (120), then SKAT0 (104); S-TWASYFS only identifies 29 genes and S-PrediXcanYFS identifies 31 genes. Meanwhile, among 2682 genome-wide significant genes for HDL2, Overall identifies the largest number of genes (192); OT and SKATO0 identify 190 genes; S-TWASMETSIM and S-PrediXcanMETSIM identify 112 and 118 genes. respectively.

Discussions

In this paper, we develop a powerful and computationally efficient method, Overall, for gene-based association studies using GWAS summary data. Overall aggregates information from three traditional types of gene-based association tests (BT, SKAT, SKATO) and also incorporates eQTL data. Both our simulation studies and real data analysis confirm that our proposed method can control type I error rates correctly and has very good performance compared with other comparison methods. In “Real data analysis”, Overall identify more significant genes than other methods, and there are some genes reported by GWAS catalog which are only identified by Overall.

There are some advantages of our proposed method. First, Overall adaptively aggregates information from multiple gene-based association tests. Most combination tests (i.e., Fisher’s combination test[61]) assume that the p-values should be calculated from independent tests. To combine information from highly correlated gene-based association tests, Overall utilizes the extended Simes procedure[5,22]. It is shown that this procedure to combine multiple tests is stable and effective regardless of whether the tests are highly correlated[24,62]. Second, Overall is more powerful than the traditional gene-based association tests, some popular transcriptome association tests (i.e., S-PrediXcan[29] and S-TWAS[12]), and other eQTL weighted combination tests (i.e., ominous test[15]). By aggregating information from different tests and incorporating multiple eQTL—derived weights, Overall can achieve a higher statistical power under a variety of situation settings. Meanwhile, our simulation studies and real data analyses show that the extended Simes procedure is more powerful than the Cauchy combination method, especially if the proportion of causal variants in a gene is small. Third, the p-values of Overall can be analytically computed without using permutations, therefore, Overall is computationally efficient. Finally, using the regularization procedure to correct the estimated LD can reduce the potential statistical noise in the LD estimation if LD is estimated using a reference panel with small sample size. In addition, Overall can be easily applied to genetic association studies with either individual-level data or GWAS summary statistics.

In this paper, we combine three types of traditional gene-based association tests (BT, SKAT, SKATO). However, the combination procedure used in the paper is very general. Other more powerful gene-based association tests can also be combined using the same approach, such as some state-of-the-art methods (i.e., S-TWAS[12], E-MAGMA[63], and SMR[64]).

In this current study, we utilize the weights derived from four single tissue gene expression studies (CMC, METSIM, NTR, YFS). Although the extended Simes procedure in Overall allows us to employ more eQTL—derived weights from a number of studies (i.e., GTEx gene expression version 8[65] et al.), there is a possibility that the noise can be increased with the increment in the number of unrelated studies. Therefore, the power of the combination tests (i.e., Overall and OT) might be attenuated. Thus, to obtain the most robust identification of phenotypic associated genes in a real data analysis with the Overall method, we suggest incorporating eQTL datasets from the most relevant tissues to the phenotype. The last but the most important thing is that population stratification can be confounded association results[66,67]. Systematic minor allele frequency difference between transcriptomic studies of different cohorts and no matching between the estimated LD structure of Genomes Project with that in the study may increase the chances of false positive findings. Therefore, we need to eliminate false positive findings possibly caused by population stratification[68,69]. When applying the Overall method, the population of GWAS summary dataset, external reference panel (i.e., 1000 Genomes Project) used to estimate LD structure, and eQTL—derived weights should be consistent.

In this study, the computational time of the proposed method is acceptable even if the estimated correlation matrix of multiple tests is obtained by the replication procedure. Meanwhile, the estimation procedure is independent of gene-based association tests, therefore we only need to perform this procedure once for each GWAS summary dataset. For example, there are a total of 29,008 gene in the 1000 Genomes Project and we use 1000 replicates to estimate the correlation matrix of multiple tests for each gene. We perform this using the high-performance computing (HPC) cluster (Intel Xeon E5—2670 2.6 GHz, 16 GB RAM). The computational time for all genes is about 36 h CPU time with 500 nodes. Then, the p-value of the proposed method can be computed analytically which is independently performed in each GWAS summary dataset. The computational time for each GWAS dataset is about 1 h CPU time with 10 nodes.

Supplementary Information.