Integration of genome-wide association studies (GWAS) and expression quantitative trait loci (eQTL) studies is needed to improve our understanding of the biological mechanisms underlying GWAS hits, and our ability to identify therapeutic targets. Gene-level association methods such as PrediXcan can prioritize candidate targets. However, limited eQTL sample sizes and absence of relevant developmental and disease context restrict our ability to detect associations. Here we propose an efficient statistical method (MultiXcan) that leverages the substantial sharing of eQTLs across tissues and contexts to improve our ability to identify potential target genes. MultiXcan integrates evidence across multiple panels using multivariate regression, which naturally takes into account the correlation structure. We apply our method to simulated and real traits from the UK Biobank and show that, in realistic settings, we can detect a larger set of significantly associated genes than using each panel separately. To improve applicability, we developed a summary result-based extension called S-MultiXcan, which we show yields highly concordant results with the individual level version when LD is well matched. Our multivariate model-based approach allowed us to use the individual level results as a gold standard to calibrate and develop a robust implementation of the summary-based extension. Results from our analysis as well as software and necessary resources to apply our method are publicly available.
Recent technological advances allow interrogation of the genome to a high level of coverage and precision, enabling experimental studies that query the effect of genotype on both complex and molecular traits. Among these, GWAS have successfully associated genetic loci to human complex traits. GWAS meta-analyses with ever increasing sample sizes allow the detection of associated variants with smaller effect sizes [1–3]. However, understanding the mechanism underlying these associations remains a challenging problem. Another approach is the study of expression quantitative trait loci (eQTLs), measuring association between genotype and gene expression. These studies provide a wealth of biological information but tend to have smaller sample sizes. A similar observation applies to QTL studies of other traits such methylation, metabolites, or protein levels. The importance of gene expression regulation in complex traits [4–7] has motivated the development of methods to integrate eQTL studies and GWAS. To examine these mechanisms we developed PrediXcan [8], which tests the mediating role of gene expression variation in complex traits. Briefly, PrediXcan tests the hypothesis that genetic variants affect phenotypes through the regulation of gene expression traits. To do that, it correlates genetically predicted gene expression and the phenotype with the idea that causal genes are likely to show a significant association. Linear prediction models of expression using genetic variation in the vicinity of the gene are trained in reference transcriptome datasets such as Genotype-Tissue Expression project (GTEx) [9]. Due to sharing of eQTLs across multiple tissues, we have shown the benefits of an agnostic scanning across all available tissues [10]. Despite the increased multiple testing burden (for Bonferroni correction, the total number of gene-tissue pairs must be used when determining the threshold), we gain considerably in number of significant genes. However, given the substantial correlation between different tissues [9], Bonferroni correction can be too stringent increasing the false negative rate. In order to aggregate evidence more efficiently, we present here a method termed MultiXcan, which tests the joint effects of gene expression variation from different tissues. Furthermore, we develop and implement a method that only needs summary statistics from a GWAS: Summary-MultiXcan (S-MultiXcan). We make our implementation publicly available to the research community in https: //github. com/hakyimlab/MetaXcan. We apply this method to simulated and real data (222 traits from the UK Biobank study [11] and 109 public GWAS) to show the performance and proper calibration of p-values. We make all of the results publicly available at https: //doi. org/10. 5281/zenodo. 1402225.
To integrate information across tissues, MultiXcan regresses the phenotype of interest on the predicted expression of a gene in multiple tissues as follows: y = μ + t 1 g 1 + t 2 g 2 + ⋯ + t p g p + e (1) where y is the n-dimensional phenotype vector, μ is an intercept term, ti is standardized predicted expression of the gene in tissue i, gi is its effect size, and e an error term with variance σ e 2; p is the number of available tissue models. We use an F-test to assess the joint significance of the regression. Expression predictions across tissues can be highly correlated. We predicted expression for individuals from the UK Biobank cohort using models trained on 44 GTEx tissues (as presented in [10]), and found a median pair-wise correlation of rp50 = 0. 56 (IQR = 0. 69) between different tissue models in a given gene, across genes (see Methods for details). To avoid numerical issues caused by collinearity, we use principal components of the predicted expression data matrix as explanatory variables, and discard the axes of smallest variation (PCA regularization). Additional covariates can be added to the regression seamlessly. Fig 1-a displays an overview of the method; see further details in the Methods section. S1 Fig shows an example of the correlation between tissues of predicted expression of the gene SLC5A6. We applied MultiXcan to 222 traits from the UK Biobank cohort. The traits were chosen based on several criteria, such as availability of well-established literature, binary traits having enough cases, or potential interest for a phenome-wide study (allergy, behavioral, metabolic and anthropometric phenotypes). We used Elastic Net prediction models trained on 44 tissues from GTEx, originally presented in [10]. We compared three approaches for assessing the significance of a gene jointly across all tissues: 1) running PrediXcan using the most relevant tissue; 2) running PrediXcan using all tissues, one tissue at a time; 3) running MultiXcan. Fig 1-b illustrates the results from each approach. We summarize a comparison between approaches 2) and 3) in Table 1. PrediXcan overcomes MultiXcan only in 21 traits, all of them with less than 50 significant associations across both methods. MultiXcan detects more associations in 103 traits. Fig 2-a and 2-b show a comparison of detections for both MultiXcan and PrediXcan. See S1 Dataset for a summary of detections per trait, and S2 and S3 Datasets for the full list of significant MultiXcan and PrediXcan results respectively. As an illustrative example, we examined more closely the results for self-reported high cholesterol phenotype (http: //biobank. ctsu. ox. ac. uk/crystal/field. cgi? id=20002). We used 50,497 cases and 100,994 controls. After Bonferroni correction, MultiXcan was able to detect a larger number of significantly associated genes (251 detections) than PrediXcan using all tissues (196 detections) or only a single tissue (whole blood, 33 detections). 172 genes were detected by both PrediXcan and MultiXcan. Fig 2-c shows the QQ-plot for associations in these three approaches. There are 79 genes associated to high cholesterol via MultiXcan and not PrediXcan. Among them, we find genes related to lipid metabolism (APOM [12], PAFAH1B2 [13]), glucose transport (SLC5A6 [14]), and vascular processes (NOTCH4 [15]). The well known gene SORT1 is detected by both MultiXcan and PrediXcan. To evaluate MultiXcan’s performance in different known scenarios, we simulated traits as a function of different numbers of causal tissues for each gene: a single tissue, multiple tissues, all available tissues. We executed PrediXcan, MultiXcan without PCA regularization, and MultiXcan with PCA regularization. We show proper calibration under the null hypothesis of no association in S3 Fig, and robustness of the regularization approach in S6 Fig. See further details in S1 Supplementary Note. As expected, when there is a known single causal tissue, PrediXcan with the known tissue yields more significant associations. However, when there are multiple causal tissues, MultiXcan yields more significant associations than the best single tissue PrediXcan results. In traits simulated from a single causal tissue, PrediXcan outperforms MultiXcan in 99. 9% of the cases (AOV p-value < 10−16). MultiXcan performs best in scenarios with multiple causal tissues (84. 4% of the times when a few tissues are causal, and 99. 5% when all tissues are causal; AOV p-value < 10−16 in both cases). One caveat is that the simulation does not cover cases when the prediction in the single tissue has low quality. In such an scenario, borrowing information from other tissues will still be beneficial. To expand the applicability of our method to massive sample sizes and to studies where individual level data are not available, we extend our method to use summary results rather than individual-level data. We call this extension Summary-MultiXcan (S-MultiXcan). We infer the joint estimates of effect sizes of predicted expression on phenotype (Eq 1) using the marginal estimates. We also compute the covariance matrix of the effect sizes and leverage the asymptotic multivariate normality of the estimates, to compute a statistic that is approximately χ p 2 (p number of tissues). The final expression is equivalent to the omnibus test mentioned in [16], which can be interpreted as a specific case of general weighted association analysis [17]. Fig 3-a illustrates our approach and the details can be found in the Methods section. As with the individual level approach, the correlation between tissues leads to numerical problems (due to near singular covariance matrices that need to be inverted). We address this by using a pseudo inverse approach which, in a nutshell, uses singular value decomposition (SVD) of the covariance matrix to keep only the components of large variation. This is analogous to the PCA regularization used for the individual level approach. Thus we test for significance using χ k 2 with k the number of surviving components. See details in the Methods Section. A robust implementation for calculating predicted expression correlation is critical to avoid unnecessary false positive results. In principle, it is possible to simply calculate the correlation between tissues using predicted expression in a reference set. However, we found that this approach can lead to large differences between the individual level data results (our gold standard) and the summary level ones when SNPs from the reference set are missing in the GWAS results. An example of this is shown in S8 Fig with the Type 1 Diabetes study from the Wellcome Trust Case-Control Consortium (WTCCC); association data is included in S9 Dataset. To avoid this problem, we calculate the covariance matrix between tissues using only the predictor SNPs that are common in both the GWAS summary and the reference LD set. Fig 4 displays a few examples of the general agreement between the individual-level MultiXcan and S-MultiXcan. The summary-based version’s results tend to be slightly more conservative than MultiXcan, as illustrated in S2 Fig. As a general comparison to the individual-level method, we list a summary of S-MultiXcan’s application to the 222 UK Biobank traits on Table 2; we observe an adequate similarity between S-MultiXcan’s and MultiXcan’s summaries. The small loss in power arises from the imperfect match of LD between the UK cohort and the reference panel. To reduce false positives due to LD misspecification when dealing with GWAS summary statistics, we discard any significant association result for a gene if the best single tissue result has p-value greater than 10−4 (“suspicious associations”). In other words, we keep significant associations if at least one single gene-tissue pair association is borderline significant or better (10−5 is the Bonferroni threshold for a typical tissue model). This is rather conservative since it is possible that evidence with modest significance from weakly correlated tissues can lead to very significant combined association when their effects get aggregated. For example among Bonferroni significant genes in the individual level analysis, a median of 8. 3% across traits (IQR = 5. 7%) have the most significant marginal (PrediXcan) p-value greater than 10−4. We list the number of such genes for each of the 222 UK Biobank traits in S8 Dataset. We applied S-MultiXcan to 109 traits on publicly available GWAS, chosen with a similar criteria as UK Biobank’s traits. Like the individual level method, we observed S-MultiXcan to detect more associations than S-PrediXcan in most cases (average detection increase 10), as shown in Fig 3-b, after discarding suspicious associations. We also show the QQ-plots for a sample trait (Schizophrenia) on Fig 3-c and the total number of associations across all public GWAS traits in 3-d. We display a summarized comparison between S-MultiXcan and S-PrediXcan in S1 Table, after discarding suspicious associations. The list of analyzed traits can be found in S4 and S5 Datasets contains a summary of significant associations for each trait and for each method. S6 Dataset lists the significant S-MultiXcan results for each trait. These results have been uploaded to https: //doi. org/10. 5281/zenodo. 1402225.
We make our software publicly available on a GitHub repository: https: //github. com/hakyimlab/MetaXcan. Prediction model weights and covariances for different tissues can be downloaded from http: //predictdb. org/. A short working example can be found on the GitHub page; more extensive documentation can be found on the project’s https: //github. com/hakyimlab/MetaXcan/wiki. The results of S-MultiXcan applied to the 44 human tissues and a broad set of phenotypes can be queried on http: //gene2pheno. org. The data used in this paper is publicly available in https: //doi. org/10. 5281/zenodo. 1402225.
This study uses de-identified genotype and phenotype data from public repositories including dbGaP, EGA, and UK Biobank. Our study has been determined to be non-human subject research by the University of Chicago’s IRB protocol number IRB16-0921. MultiXcan consists of fitting a linear regression of the phenotype on predicted expression from multiple tissue models jointly: y = ∑ j = 1 p t j g j + e = T g + e (2) where y is a centered vector of phenotypes for n individuals, tj is an n-vector of standardized predicted gene expression for model j, gj is the effect size for the predicted gene expression j, e is an error term with variance σ e 2, and p is the number of tissues. Thus, T is a data matrix where each column j contains the values from tj, and g is the p-vector of effect sizes gj. The high degree of eQTL sharing between different tissues induces a high correlation between predicted expression levels. In order to avoid collinearity issues and numerical instability, we decompose the predicted expression matrix into principal components and keep only the eigenvectors of non negligible variance. To select the number of components, we used a condition number threshold of λ max λ i < 30, where λi is an eigenvalue of the matrix Tt T. As a side effect, we observe moderate increases in significance levels because less informative components of tissue expression are discarded from the model. A range of values between 10 and 100 yielded similar results in the simulations described in S1 Supplementary Note as displayed in S6 Fig. Lastly, we use an F-test to quantify the significance of the joint fit. We use Bonferroni correction to determine the significance threshold. For MultiXcan, we use the total number of genes with a prediction model in at least one tissue, which yields a threshold approximately at 0. 05/17500 ∼ 2. 9 × 10−6. For PrediXcan across all tissues, we use the total number of gene-tissue pairs, which yields a threshold approximately at 0. 05/200,000 ∼ 2. 5 × 10−7. Since the tested hypotheses are not independent, Bonferroni correction is overly conservative, as can be seen when counting the number of associations via FDR in S7 Fig. We have demonstrated that S-PrediXcan can accurately infer PrediXcan results from GWAS Summary Statistics and LD information from a reference panel [10], with the added benefits of reduced computational and regulatory burden. Here we extend MultiXcan in a similar fashion. Summary-MultiXcan (S-MultiXcan) infers the individual-level MultiXcan results, using univariate S-PrediXcan results and LD information from a reference panel. It consists of the following steps: Prediction Models were obtained from http: //predictdb. org/ resource. These models were trained using Elastic Net as implemented in R’s package glmnet [38], with a mixing parameter α = 0. 5, on 44 tissue studies from GTEx’ release version 6p. The underlying GTEx study data was obtained from dbGaP with accesion number phs000424. v6. p1. Please see [10] for details. We implemented MultiXcan and S-MultiXcan using python scientific packages, working up from existing software in the MetaXcan package. S-PrediXcan, PrediXcan, MultiXcan and S-MultiXcan analysis were computed using the Center for Research Informatics’ high performance cluster at the University of Chicago. PrediXcan, S-PrediXcan, MultiXcan and S-MultiXcan results have been uploaded to the http: //gene2pheno. org resources. The databases are open to the research community for arbitrary programmatic query.
We develop a new method, MultiXcan, to test the mediating role of gene expression variation on complex traits, integrating information available across multiple tissue studies. We show this approach has higher power than traditional single-tissue methods. We extend this method to use only summary-statistics from public GWAS. We apply these methods to 222 complex traits available in the UK Biobank cohort, and 109 complex traits from public GWAS and discuss the findings.
genome-wide association studies
statistics
sociology
computational biology
social sciences
random variables
covariance
singular value decomposition
mathematics
forecasting
algebra
genome analysis
research and analysis methods
mathematical and statistical techniques
gene expression
consortia
probability theory
linear algebra
phenotypes
genetics
biology and life sciences
physical sciences
genomics
gene prediction
statistical methods
human genetics
Integrating predicted transcriptome from multiple tissues improves association detection