Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics

1Over the last decade, GWAS have been successful in robustly associating genetic loci to human complex traits. However, the mechanistic understanding of these discoveries is still limited, hampering the translation of the associations into actionable targets. Studies of enrichment of expression quantitative trait loci (eQTLs) among trait-associated variants– show the importance of gene expression regulation. Functional class quantification showed that 80% of the common variant contribution to phenotype variability in 12 diseases can be attributed to DNAase I hypersensitivity sites, further highlighting the importance of transcript regulation in determining phenotypes.

2Many transcriptome studies have been conducted where genotypes and expression levels are assayed for a large number of individuals–. The most comprehensive transcriptome dataset, in terms of examined tissues, is the Genotype-Tissue Expression Project (GTEx): a large-scale effort where DNA and RNA were collected from multiple tissue samples from nearly 1000 individuals and sequenced to high coverage,. This remarkable resource provides a comprehensive cross-tissue survey of the functional consequences of genetic variation at the transcript level.

3To integrate knowledge generated from these large-scale transcriptome studies and shed light on disease biology, we developed PrediXcan, a gene-level association approach that tests the mediating effects of gene expression levels on phenotypes. PrediXcan is implemented on GWAS or sequencing studies (i.e., studies with genome-wide interrogation of DNA variation and phenotypes). It imputes transcriptome levels with models trained in measured transcriptome datasets (e.g., GTEx). These predicted expression levels are then correlated with the phenotype in a gene association test that addresses some of the key limitations of GWAS.

4Meta-analysis efforts that aggregate results from multiple GWAS have been able to identify an increasing number of associations that were not detected with smaller sample sizes–. We will refer to these results as Genome-wide association meta-analysis (GWAMA) results. In order to harness the power of these increased sample sizes while keeping the computational burden manageable, methods that use summary level data rather than individual level data are needed.

5Methods similar to PrediXcan that estimate the association between intermediate gene expression levels and phenotypes, but use summary statistics have been reported: TWAS (summary version) and Summary Mendelian Randomization (SMR). Another class of methods that integrate eQTL information with GWAS results are based on colocalization of eQTL and GWAS signals. Colocalized signals provide evidence of possible causal relationship between the target gene of an eQTL and the complex trait. These include RTC, Sherlock, COLOC, and more recently eCAVIAR and ENLOC.

6Here we derive a mathematical expression that allows us to compute the results of PrediXcan without the need to use individual-level data, greatly expanding its applicability. We compare with existing methods and outline a best practices framework to perform integrative gene mapping studies, which we term MetaXcan.

7We apply the MetaXcan framework by first training over one million elastic net prediction models of gene expression traits, covering protein coding genes across 44 human tissues from GTEx, and then performing gene-level association tests over 100 phenotypes from 40 large meta-analysis consortia and dbGaP.

8We have derived an analytic expression to compute the outcome of PrediXcan using only summary statistics from genetic association studies. Details of the derivation are shown in the Methods section. In Fig. 1a we illustrate the mechanics of Summary-PrediXcan (S-PrediXcan) in relation to traditional GWAS and the individual-level PrediXcan method.

9We find high concordance between PrediXcan and S-PrediXcan results indicating that in most cases, we can use the summary version without loss of power to detect associations. Figure 2 shows the comparison of PrediXcan and S-PrediXcan Z-scores for a simulated phenotype (under the null hypothesis), a cellular growth phenotype and two disease phenotypes: type 1 diabetes and bipolar disorder from the WTCCC Consortium; see Supplementary Notes 1, 2 and 3 for details. For the simulated phenotype, the study sets (in which GWAS is performed) and the reference set (in which LD between SNPs is computed) were African, East Asian, and European subsets from 1000 Genomes. The training set (in which prediction models are trained) was European (DGN Cohort) in all cases. The high correlation between PrediXcan and S-PrediXcan demonstrates the robustness of our method to mismatches between reference and study sets. Despite the generally good concordance between the summary and individual level methods, there were a handful of false positive results with S-PrediXcan much more significant than PrediXcan. This underscores the need to use closely matched LD information whenever possible. Supplementary Fig. 11 shows S-PrediXcan’s performance on a phenotype simulated under the alternative hypothesis.

10Notice that we are not testing here whether PrediXcan itself is robust to population differences between training and study sets. Robustness of the prediction across populations has been previously reported. We further corroborated this in Supplementary Fig. 10.

11Next we compare with other summary result-based methods such as S-TWAS, SMR, and COLOC.

12One class of methods seeks to determine whether eQTL and GWAS signals are colocalized or are distinct although linked by LD. This class includes COLOC, Sherlock, and RTC, and more recently eCAVIAR, and ENLOC. Thorough comparison between these methods can be found in refs.,. HEIDI, the post filtering step in SMR that estimates heterogeneity of GWAS and eQTL signals, can be included in this class. We focus here on COLOC, whose quantification of the probability of five configurations complements well with S-PrediXcan results.

13COLOC provides the probability of five hypotheses: H0 corresponds to no eQTL and no GWAS association, H1 and H2 correspond to association with eQTL but no GWAS association or vice-versa, H3 corresponds to eQTL and GWAS association but independent signals, and finally H4 corresponds to shared eQTL and GWAS association. P0, P1, P2, P3, and P4 are the corresponding probabilities for each configuration. The sum of the five probabilities is 1. The authors recommend to interpret H0, H1, and H2 as limited power; for convenience we will aggregate these three hypotheses into one event with probability 1-P3-P4.

14Figure 3 shows ternary plots with P3, P4, and 1-P3-P4 as vertices. The blue region, top subtriangle, corresponds to high probability of colocalized eQTL and GWAS signals (P4). The orange region at bottom left corresponds to high probability of distinct eQTL and GWAS signals (P3). The gray region at center and bottom right corresponds to low probability of both colocalization and independent signals.

15Figure 3b shows association results for all gene-tissue pairs with the height phenotype. We find that most associations fall in the gray, “undetermined,” region. When we restrict the plot to S-PrediXcan Bonferroni-significant genes (Fig. 3c), three distinct peaks emerge in the high P4 region (P4 > 0.5, “colocalized signals”), high P3 region (P3 > 0.5, “independent signals” or “non-colocalized signals”), and “undetermined” region. Moreover, when genes with low prediction performance are excluded (Supplementary Fig. 6d) the “undetermined” peak significantly diminishes.

16These clusters provide a natural way to classify significant genes and complement S-PrediXcan results. Depending on false positive/false negative trade-off choices, genes in the “independent signals” or both “independent signals” and “undetermined” can be filtered out. The proportion of colocalized associations (P4 > 0.5) ranged from 5 to 100% depending on phenotype with a median of 27.6%. The proportion of “non-colocalized” associations ranged from 0 to 77% with a median of 27.0%. Supplementary Table 1 summarizes the percentages of significant associations that fall into the different colocalization regions.

17This post-filtering idea was first implemented in the SMR approach using HEIDI. Comparison of COLOC results with HEIDI is shown in Supplementary Fig. 6e to 6h.

18Gusev et al. have proposed Transcriptome-Wide Association Study based on summary statistics (STWAS), which imputes the SNP level Z-scores into gene level Z-scores. This is not the same as computing the results of individual level TWAS. We show (in Methods section) that the difference between the individual level and summary level TWAS is given by the factor, where Rl is the proportion of variance in the phenotype explained by a SNP’s allelic dosage, and Rg is the proportion explained by gene expression (see Methods section). For most practical purposes we have found that this factor is very close to 1 so that if the same prediction models were used, no substantial difference between S-TWAS and S-PrediXcan should be expected.

19Figure 4a shows a diagram of S-PrediXcan and S-TWAS. Both use SNP to phenotype associations results (ZX,Y) and prediction weights (wX,Tg) to infer the association between the gene expression level (Tg) and phenotype (Y).

20Figure 4b compares S-TWAS significance (as reported in ref. ) to S-PrediXcan significance. The difference between the two approaches is mostly driven by the different prediction models: TWAS uses BSLMM whereas PrediXcan uses elastic net. BSLMM allows two components: one sparse (small set of large effect predictors) and one polygenic (all variants contribute some marginal effect to the prediction). For PrediXcan we have chosen to use a sparse model (elastic net) based on the finding that the genetic component of gene expression levels is mostly sparse.

21Figure 4c shows that the COLOC-estimated proportion of non-colocalized (independent) GWAS and eQTL signals is larger among TWAS significant genes than among S-PrediXcan significant ones. We believe this is due to the polygenic component of BSLMM models, a wider set of SNPs increasing the chances of COLOC yielding a non-colocalized result. Figure 4d shows that, for most traits, the COLOC-estimated proportion of colocalized signals is larger among S-PrediXcan significant genes than S-TWAS significant genes.

22Zhu et al. have proposed Summary Mendelian Randomization (SMR), a summary data based Mendelian randomization that integrates eQTL results to determine target genes of complex trait-associated GWAS loci. They derive an approximate -statistic (Eq 5 in ref. ) for the mediating effect of the target gene expression on the phenotype. Figure 5a depicts this mechanism.

23Unfortunately, the derived statistic is not well calibrated. A QQ plot comparing the SMR statistic (under the null hypothesis of genome-wide significant eQTL signal and no GWAS association) shows deflation. The sample mean of the statistic is ≈0.93 instead of 1, the expected value for the mean of a random variable (see Fig. 5e, f and Methods section for details). The approximation is only valid in two extreme cases: when the eQTL association is much stronger than the GWAS association or vice versa, when the GWAS association is much stronger than the eQTL association (see Methods section for details).

24One limitation is that the significance of the SMR statistic is the lower of the top eQTL association (genotype to expression) or the GWAS association (genotype to phenotype) as shown in Fig. 5e, f. Given the much larger sample sizes of GWAS studies, for most genes, the combined significance will be determined by the eQTL association. The combined statistic forces us to apply multiple testing correction for all genes, even those that are distant to GWAS associated loci, which is unnecessarily conservative. Keep in mind that currently both SMR and PrediXcan only use cis associations. An example may clarify this further. Let us suppose that for a given phenotype there is only one causal SNP and that the GWAS yielded a highly significant p-value, say 10−50. Let us also suppose that there is only one gene (gene A) in the vicinity (we are only using cis predictors) associated with the causal SNP with p = 10−5. SMR would compute the p-values of all genes and yield a p-value ≈ 10−5 for gene A (the less significant p-value). However, after multiple correction this gene would not be significantly associated with the phenotype. Here it is clear that we should not be adjusting for testing of all genes when we know a priori that only one is likely to produce a gene level association. In contrast, the PrediXcan p-value would be ≈10−50 for gene A and would be distributed uniformly from 0 to 1 for the remaining genes. Most likely only gene A (or perhaps a handful of genes, just by chance) would be significant after Bonferroni correction. If we further correct for prediction uncertainty (here = eQTL association), a p-value of ≈10−5 would remain significant since we only need to correct for the (at most) handful of genes that were Bonferroni significant for the PrediXcan p-value.

25Another potential disadvantage of this method is that only top-eQTLs are used for testing the gene level association. This does not allow to aggregate the effect on the gene across multiple variants.

26Figure 5b compares S-PrediXcan (elastic net) and SMR association results. As expected, SMR p-values tend to be less significant than S-PrediXcan’s in large part due to the additional adjustment for the uncertainty in the eQTL association. Figure 5c, d show that the SMR significance is bounded by the eQTL and GWAS association strengths of the top eQTL. Figure 5g shows a comparison between SMR’s and S-PrediXcan’s proportion of non-colocalization, while Fig. 5h compares proportion of colocalization, as estimated by COLOC. SMR shows a higher proportion of colocalized and independent signals. This is expected since SMR uses a more stringent eQTL association criterion so that there are few significant genes in the undetermined region.

27SMR introduces a post filtering step via an approach called HEIDI, which is compared to COLOC in Fig. 3e and Supplementary Fig. 6.

28Building on S-PrediXcan and existing approaches, we define a general framework (MetaXcan) to integrate eQTL information with GWAS results and map disease-associated genes. This evolving framework can incorporate models and methods to increase the power to detect causal genes and filter out false positives. Existing methods fit within this general framework as instances or components (Fig. 6a).

29The framework starts with the training of prediction models for gene expression traits followed by a selection of high-performing models. Next, a mathematical operation is performed to compute the association between each gene and the downstream complex trait. Additional adjustment for the uncertainty in the prediction model can be added. To avoid capturing LD-contaminated associations, which can occur when expression predictor SNPs and phenotype causal SNPs are different but in LD, we use colocalization methods that estimate the probability of shared or independent signals.

30PrediXcan implementations use elastic net models motivated by our observation that gene expression variation is mostly driven by sparse components. TWAS implementations have used Bayesian Sparse Linear Mixed Models (BSLMM). SMR fits into this scheme with prediction models consisting solely of the top eQTL for each gene (weights are not necessary here since only one SNP is used at a time).

31For the last step, we chose COLOC to estimate the probability of colocalization of GWAS and eQTL signals. COLOC probabilities cluster more distinctly into different classes and thus, unlike other methods, suggests a natural cut off threshold at P = 0.5. Another advantage of COLOC is that for genes with low probability of colocalization, it further distinguishes distinct GWAS and eQTL signals from low power. This is a useful feature that future development of colocalization methods should also offer. SMR, on the other hand, uses its own estimate of “heterogeneity” of signals calculated by HEIDI.

32Perform PrediXcan or S-PrediXcan using all tissues. Use Bonferroni correction for all gene-tissue pairs: keep p < 0.05/number of gene-tissue pairs tested.Keep associations with significant prediction performance adjusting for number of PrediXcan significant gene-tissue pairs: keep prediction performance p-values < 0.05/(number of significant associations from previous step).Filter out LD-contaminated associations, i.e., gene-tissue pairs in the “independent signal” (=”non-colocalized”) region of the ternary plot (See Fig. 3a): keep COLOC P3 < 0.5 (Blue and gray regions in Fig. 3a).If further reduction of number of genes to be taken to replication or validation is desired, keep only hits with explicit evidence of colocalization: P4 > 0.5 (Blue region in Fig. 3a).

33Perform PrediXcan or S-PrediXcan using all tissues. Use Bonferroni correction for all gene-tissue pairs: keep p < 0.05/number of gene-tissue pairs tested.

34Keep associations with significant prediction performance adjusting for number of PrediXcan significant gene-tissue pairs: keep prediction performance p-values < 0.05/(number of significant associations from previous step).

35Filter out LD-contaminated associations, i.e., gene-tissue pairs in the “independent signal” (=”non-colocalized”) region of the ternary plot (See Fig. 3a): keep COLOC P3 < 0.5 (Blue and gray regions in Fig. 3a).

36If further reduction of number of genes to be taken to replication or validation is desired, keep only hits with explicit evidence of colocalization: P4 > 0.5 (Blue region in Fig. 3a).

37Any choice of thresholds has some level of arbitrariness. Depending on the false positive and negative trade off, these numbers may be changed.

38We downloaded summary statistics from meta analyses of over 100 phenotypes from 40 consortia. The full list of consortia and phenotypes is shown in Supplementary Data 2. We tested association between these phenotypes and the predicted expression levels using elastic net models in 44 human tissues from GTEx as described in the Methods section, and a whole blood model from the DGN cohort presented in ref.. We illustrate this application in Fig. 6b.

39S-PrediXcan’s results tend to be more significant as the genetic component of gene expression increases (larger cross-validated prediction performance ). Similarly, S-PrediXcan associations tend to be more significant when prediction is more reliable (p-values of association between predicted and observed expression levels are more significant, i.e., when prediction performance p-value is smaller). The trend is seen both when results are averaged across all tissues for a given phenotype or across all phenotypes for a given tissue, as displayed in Supplementary Figs. 1-4. This trend was also robust across different monotonic functions of the Z-scores.

40We used a Bonferroni threshold accounting for all the gene-tissue pairs that were tested (0.05/total number of gene-tissue pairs ≈2.5e-7). This approach is conservative because the correlation between tissues would make the total number of independent tests smaller than the total number of gene-tissue pairs. Height had the largest number of significantly associated unique genes at 1686 (based on a GWAMA of 250 K individuals). Other polygenic diseases with a large number of associations include schizophrenia with 305 unique significant genes (n = 150 K individuals), low-density lipoprotein cholesterol (LDL-C) levels with 296 unique significant genes (n = 188 K), other lipid levels, glycemic traits, and immune/inflammatory disorders such as rheumatoid arthritis and inflammatory bowel disease. For other psychiatric phenotypes, a much smaller number of significant associations was found, with eight significant genes for bipolar disorder (n = 16,731) and one for major depressive disorder (n = 18,759), probably due to smaller sample sizes, but also smaller effect sizes.

41When step 2 from the suggested pipeline is applied, keeping only reliably predicted genes, we are left with 739 genes for height, 150 for schizophrenia, 117 for LDL-C levels.

42After step 3, which keeps genes that are without strong evidence of LD-contamination, these numbers dropped to 264 for height, 58 for schizophrenia, and 60 for LDL-C levels. After step 4, which keeps only genes with strong evidence of colocalization, we find 215 genes for height, 49 for schizophrenia, and 35 for LDL-C. The counts for the full set of phenotypes can be found in Supplementary Data 2.

43Mostly, genome-wide significant genes tend to cluster around known SNP-level genome-wide significant loci or sub-genome-wide significant loci. Regions with sub-genome-wide significant SNPs can yield genome-wide significant results in S-PrediXcan, because of the reduction in multiple testing and the increase in power arising from the combined effects of multiple variants. Supplementary Table 2 lists a few examples where this occurs.

44The proportion of colocalized associations (P4 > 0.5) ranged from 5 to 100% depending on phenotype with a median of 27.6%. The proportion of “non colocalized” associations ranged from 0 to 77% with a median of 27.0%.

45See full set of results in our online catalog (gene2pheno.org). Significant gene-tissue pairs are included in Supplementary Data 3. To facilitate comparison, the catalog contains all SMR results we generated and the S-TWAS results reported by ref. for 30 GWAS traits and GTEx BSLMM models. Note that SMR application to 28 phenotypes was reported by ref. using whole blood eQTL results from ref..

46We reasoned that if complete knock out of monogenic disease genes cause severe forms of the disease, more moderate alterations of gene expression levels (as affected by regulatory variation in the population) could cause more moderate forms of the disease. Thus moderate alterations in expression levels of monogenic disease genes (such as those driven by eQTLs) may have an effect on related complex traits, and this effect could be captured by S-PrediXcan association statistics. To test this hypothesis, we obtained genes listed in the ClinVar database for obesity, rheumatoid arthritis, diabetes, Alzheimer’s, Crohn’s disease, ulcerative colitis, age-related macular degeneration, schizophrenia, and autism. Figure 7 displays the QQ plot for all associations and compares to those in ClinVar database. As postulated, we found enrichment of significant S-PrediXcan associations for ClinVar genes for all tested phenotypes except for autism and schizophrenia. The lack of significance for autism is probably due to insufficient power: the distribution of p-values is close to the null distribution. In contrast, for schizophrenia, many genes were found to be significant in the S-PrediXcan analysis. There are several reasons that may explain this lack of enrichment: genes identified with GWAS and subsequently with S-PrediXcan have rather small effect sizes, so that it would not be surprising that they were missed until very large sample sizes were aggregated; ClinVar genes may originate from rare mutations that are not well covered by our prediction models, which are based on common variation (due to limited sample sizes of eQTL studies and the minor allele frequency –MAF filter used in GWAS studies); or the mechanism of action of the schizophrenia linked ClinVar genes may be different than the alteration of expression levels. Also, the pathogenicity of some of the ClinVar entries has been questioned. The list of diseases in ClinVar used to generate the enrichment figures can be found in Supplementary Data 1, along with the corresponding association results.

47Most genes were found to be significantly associated in a handful of tissues as illustrated in Fig. 8b. For example, for LDL-C levels, liver was the most enriched tissue in significant associations as expected given known biology of this trait (See Supplementary Fig. 5). This prominent role of liver was apparent despite the smaller sample size available for building liver models (n = 97), which was less than a third of the numbers available for muscle (n = 361) or lung (n = 278).

48However, in general, tissues expected to stand out as more enriched for diseases given currently known biology did not consistently do so when we looked at the average across all (significant) genes, using various measures of enrichment. For example, the enrichment in liver was less apparent for high-density lipoprotein cholesterol (HDL-C) or triglyceride levels. We find for many significant associations that the evidence is present across multiple tissues. This may be caused by a combination of context specificity and sharing of regulatory mechanism across tissues.

49Next, we illustrate the challenges of identifying disease relevant tissues based on eQTL information using three genes with well established biology: C4A for schizophrenia and SORT1 and PCSK9 both for LDL-C and cardiovascular disease. S-PrediXcan results for these genes and traits, and regulatory activity by tissue (as measured by the proportion of expression explained by the genetic component), are shown in Fig. 8a. Representative results are shown in Supplementary Tables 3, 4 and 5. Supplementary Data 4 contains the full set of MetaXcan results (i.e., association, colocalization, and HEIDI) for these genes.

50SORT1 is a gene with strong evidence for a causal role in LDL-C levels, and as a consequence, is likely to affect risk for cardiovascular disease. This gene is most actively regulated in liver (close to 50% of the expression level of this gene is determined by the genetic component) with the most significant S-PrediXcan association in liver (p-value ≈ 0, Z = −28.8), consistent with our prior knowledge of lipid metabolism. In this example, tissue specific results suggest a causal role of SORT1 in liver.

51However, in the following example, association results across multiple tissues do not allow us to discriminate the tissue of action. C4A is a gene with strong evidence of causal effect on schizophrenia risk via excessive synaptic pruning in the brain during development. Our results show that C4A is associated with schizophrenia risk in all tissues (p < 2.5 × 10−7 in 36 tissue models and p < 0.05 for the remaining four tissue models).

52PCSK9 is a target of several LDL-C lowering drugs currently under trial to reduce cardiovascular events. The STARNET study profiled gene expression levels in cardiometabolic disease patients and showed tag SNP rs12740374 to be a strong eQTL for PCSK9 in visceral fat but not in liver. Consistent with this, our S-PrediXcan results also show a highly significant association between PCSK9 and LDL-C (p ≈ 10−13) in visceral fat and not in liver (our training algorithm did not yield a prediction model for PCSK9, i.e., there was no evidence of regulatory activity). In our results, however, the statistical evidence is much stronger in tibial nerve (p ≈ 10−27). Accordingly, in our training set (GTEx), there is much stronger evidence of regulation of this gene in tibial nerve compared to visceral fat.

53Most associations highlighted here have high colocalization probabilities. See Supplementary Tables 3, 4, and 5. However, visceral fat association shows evidence of non colocalization (probability of independent signals P3 = 0.69 in LDL-C). It is possible that the relevant regulatory activity in visceral adipose tissue was not detected in the GTEx samples for various reasons but it was detected in tibial nerve. Thus by looking into all tissues’ results we increase the window of opportunities where we can detect the association.

54PCSK9 yields colocalized signals for LDL-C levels in Tibial Nerve, Lung, and Whole blood. SORT1 shows colocalization with LDL-C in liver (P4 ≈ 1) and pancreas (P4 = 0.90). C4A is colocalized with schizophrenia risk for the majority of the tissues (29/40) with a median colocalization probability of 0.82.

55These examples demonstrate the power of studying regulation in a broad set of tissues and contexts and emphasize the challenges of determining causal tissues of complex traits based on in-silico analysis alone. Based on these results, we recommend to scan all tissues’ models to increase the chances to detect the relevant regulatory mechanism that mediates the phenotypic association. False positives can be controlled by Bonferroni correcting for the additional tests.

56We used data from the Resource for Genetic Epidemiology Research on Adult Health and Aging study (GERA, phs000674.v1.p1),. This is a study led by the Kaiser Permanente Research Program on Genes, Environment, and Health (RPGEH) and the UCSF Institute for Human Genetics with over 100,000 participants. We downloaded the data from dbGaP and performed GWAS followed by S-PrediXcan analysis of 22 conditions available in the European subset of the cohort.

57For replication, we chose Coronary Artery Disease (CAD), LDL cholesterol levels, Triglyceride levels, and schizophrenia, which had closely related phenotypes in the GERA study and had a sufficiently large number of Bonferroni significant associations in the discovery set. Analysis and replication of the type 2 diabetes phenotype can be found in ref.. Coronary artery disease hits were compared with “Any cardiac event,” LDL cholesterol and triglyceride level signals were compared with “Dyslipidemia,” and schizophrenia was compared to “Any psychiatric event” in GERA.

58High concordance between discovery and replication is shown in Fig. 9 where dyslipidemia association Z-scores are compared to LDL cholesterol Z-scores. The majority of gene-tissue pairs (92%, among the ones with Z-score magnitude greater than 2 in both sets) have concordant direction of effects in the discovery and replication sets. The high level of concordance is supportive of an omnigenic trait architecture.

59Following standard practice in meta-analysis, we consider a gene to be replicated when the following three conditions are met: the p-value in the replication set is <0.05, the direction of discovery and replication effects are the same, and the meta analyzed p-value is Bonferroni significant with the discovery threshold.

60We display summary statistics for this replication analysis in Table 1. Among the 56 genes significantly associated with CAD in the discovery set, 6 (11%) were significantly associated with “Any cardiac event” in GERA. Using “Dyslipidemia” as the closest matching phenotype, 78.5% and 43.5% of LDL and triglyceride genes replicated, respectively. Among the 285 genes associated with schizophrenia in the discovery set, 51 (21%) replicated. The low replication rate for CAD and Schizophrenia is likely due to the broad phenotype definitions in the replication.

61We found no consistent replication pattern difference between colocalized and non-colocalized genes.

62This is not unexpected if the LD pattern is similar between discovery and replication sets.

63The full list of significant genes can be queried in gene2pheno.org.

64Here we derive a mathematical expression to compute PrediXcan results without using individual level data, which greatly expands its applicability and is robust to study and reference set mismatches. This has not been done before. TWAS, which for the individual level approach only differs from PrediXcan on the prediction model used in the implementation, has been extended to use summary level data. When Gaussian imputation is used, the relationship between individual level and summary versions of TWAS is clear. This is not the case when extended to general weights (such as BSLMM). Our mathematical derivation shows the analytic difference between them explicitly.

65We also add a post filtering step, to mitigate issues with LD-contamination. Based on consistency with PrediXcan and interpretability of results, we have chosen to use COLOC for filtering. COLOC has the limitation of assuming a single causal variant, and has reduced power in the presence of multiple causal variants. However, colocalization estimation is an active area of research and improved versions or methods will be adopted in the future. We find that BSLMM-based TWAS results have a larger proportion of non-colocalized genes as estimated by COLOC. This could be due to the single variant assumption in COLOC but we believe this is rather a consequence of the polygenic component of BSLMM predictors. Given the predominantly sparse architecture of gene expression traits, we believe that adding a polygenic component unnecessarily increases the exposure to LD-contamination.

66Despite the generally good concordance between the summary and individual level methods, there were a handful of false positive results with S-PrediXcan much more significant than PrediXcan. This underscores the need to use closely matched LD information whenever possible.

67We applied our framework to over 100 phenotypes using transcriptome prediction models trained in 44 tissues from the GTEx Consortium and generated a catalog of downstream phenotypic association results of gene expression variation, a growing resource for the community.

68The enrichment of monogenic disease genes among related phenotype associations suggests that moderate alteration of expression levels as affected by common genetic variation may cause a continuum of phenotypic changes. Alternatively, a more complex interplay between common and rare variation could be taking place such as higher tolerance to loss of function mutations in lower expressing haplotypes which could induce association with predicted expression.

69We are finding that most trait associations are tissue specific; i.e., they are detected in a handful of tissues. However, we also find that expected tissues given known biology do not necessarily rank among the top enriched tissues. This suggests context specificity of the pathogenic mechanism; specific developmental stage or environmental conditions may be necessary to detect the regulatory event. On the other hand, we are detecting associations in unexpected tissues which suggests a sharing of regulation across multiple tissues/contexts or perhaps novel biology that takes place in these tissues. In either case, agnostic scanning of a broad set of tissues is necessary to discover these mechanisms.

70Figure 1b shows the main analytic expression used by Summary-PrediXcan for the Z-score (Wald statistic) of the association between predicted gene expression and a phenotype. The input variables are the weights used to predict the expression of a given gene, the variance and covariances of the markers included in the prediction, and the GWAS coefficient for each marker. The last factor in the formula can be computed exactly in principle, but we would need additional information that is unavailable in typical GWAS summary statistics output such as phenotype variance and sample size. Dropping this factor from the formula does not affect the accuracy of the results as demonstrated in the close to perfect concordance between PrediXcan and Summary-PrediXcan results on the diagonal of Fig. 2a.

71The approximate formula we use is:where wlg is the weight of SNP l in the prediction of the expression of gene g; is the GWAS regression coefficients for SNP l; se is standard error of, is the estimated variance of SNP l, and is the estimated variance of the predicted expression of gene g; and dosage and alternate allele are assumed to be the same.

72The inputs are based, in general, on data from three different sources: study set (e.g., GWAS study set), expression training set (e.g., GTEx, DGN), population reference set (e.g., the training set or 1000 Genomes).

73The study set is the main dataset of interest from which the genotype and phenotypes of interest are gathered. The regression coefficients and standard errors are computed based on individual-level data from the study set or a SNP-level meta-analysis of multiple GWAS. Training sets are the reference transcriptome datasets used for the training of the prediction models (GTEx, DGN, Framingham, etc.) thus the weights wlg are computed from this set. Training sets can also be used to generate variance and covariances of genetic markers, which will usually be different from the study sets. When individual level data are not available from the training set we use population reference sets such as 1000 Genomes data. In the most common use scenario, users will need to provide only GWAS results using their study set. The remaining parameters are pre-computed, and published in PredictDB.

74We display the enrichment for selected phenotypes in Supplementary Fig. 5, measured as mean(Z2). For visualization purposes, we selected 25 phenotypes from different categories such as anthropometric traits, cardiometabolic traits, autoimmune diseases, and psychiatric conditions (please see figure caption for the list of selected phenotypes). The simple mean of Z2 for all gene-tissue pairs in a phenotype was taken.

75The goal of summary-PrediXcan is to infer the results of PrediXcan using only GWAS summary statistics. Individual level data are not needed for this algorithm. We will introduce some notations for the derivation of the analytic expressions of S-PrediXcan.

76Y is the n-dimensional vector of phenotype for individuals i = 1, n. Xl is the allelic dosage for SNP l. Tg is the predicted expression (or estimated GREx, genetically regulated expression). wlg are weights to predict expression, derived from an independent training set.

77We model the phenotype as linear functions of Xl and Tgwhere α1 and α2 are intercepts, η and error terms independent of Xl and Tg, respectively. Let and be the estimated regression coefficients of Y regressed on Tg and Xl, respectively. is the result (effect size for gene g) we get from PrediXcan whereas is the result from a GWAS for SNP l.

78We will denote as and the operators that compute the sample variance and covariance, i.e.,: with. Let, and, where X′ is the p × n matrix of SNP data and is a n × p matrix where column l has the column mean of Xl (p being the number of SNPs in the model for gene g, typically p ≪ n).

79With this notation, our goal is to infer PrediXcan results ( and its standard error) using only GWAS results ( and their standard error), estimated variances of SNPs, estimated covariances between SNPs in each gene model (Γg), and prediction model weights wlg.

80Input:, se,,, wlg. Output: /se.

81Next we list the properties and definitions used in the derivationand

82The proportion of variance explained by the covariate (Tg or Xl) can be expressed as

83By definition

84Thus can be computed aswhere Wg is the vector of for SNPs in the model of g. By definition, is, the sample covariance of Xg, so that we arrive to

85Calculation of regression coefficient

86can be expressed aswhere we used the linearity of in the last step. Using Eq. (3), we arrive to

87Calculation of standard error of

88Also from the properties of linear regression we know that

89In this equation, is not necessarily known but can be estimated using the equation analogous to (6) for βl

90Thus:

91Notice that the right hand side of (8) is dependent on the SNP l while the left hand side is not. This equality will hold only approximately in our implementation since we will be using approximate values for, i.e., from reference population, not the actual study population.

92To assess the significance of the association, we need to compute the ratio of the estimated effect size and standard error se, or Z-score,

93with which we can compute the p-value as p = 2Φ(−|Zg|) where Φ(.) is the normal CDF function. Thuswhere we used Eqs. (5) and (6) in the second line and Eq. (8) in the last step. So

94Based on results with actual and simulated data for realistic effect size ranges, we have found that the last approximation does not affect our ability to identify the association. The approximation becomes inaccurate only when the effect sizes are very large. But in these cases, the small decrease in statistical efficiency induced by the approximation is compensated by the large power to detect the larger effect sizes.

95The variance of predicted expression is computed using Eq. (4) which takes weights for each SNP in the prediction model and the correlation (LD) between the SNPs. The correlation is computed in a reference set such as 1000G or in the training set.

96To train our prediction models, we obtained genotype data and normalized gene expression data collected by the GTEx Project. We used 44 different tissues sampled by GTEx and thus generated 44 different tissue-wide models (dbGaP Accession phs000424.v6.p1). Sample sizes for different tissues range from 70 (Uterus) to 361 (Muscle—Skeletal). The models referenced in this paper make use of the GTEx Project’s V6p data, a patch to the version 6 data and makes use of improved gene-level annotation. We removed ambiguously stranded SNPs from genotype data, i.e., ref/alt pairs A/T, C/G, T/A, G/C. Genotype data was filtered to include only SNPs with MAF > 0.01. For each tissue, normalized gene expression data was adjusted for covariates such as gender, sequencing platform, the top three principal components from genotype data and top PEER Factors. The number of PEER Factors used was determined by sample size: 15 for n < 150, 30 for n between 150 and 250, and 35 for n > 250. Covariate data was provided by GTEx. For our analysis, we used protein-coding genes only.

97For each gene-tissue pair for which we had adjusted expression data, we fit an Elastic-Net model based on the genotypes of the samples for the SNPs located within 1 Mb upstream of the gene’s transcription start site and 1 Mb downstream of the transcription end site. We used the R package glmnet with mixing parameter alpha equal to 0.5, and the penalty parameter lambda was chosen through 10-fold

98Once we fit all models, we retained only those with q-value less than 0.05 For each tissue examined, we created a sqlite database to store the weights of the prediction models, as well as other statistics regarding model training. Supplementary Table 6 contains summary statistics on the models for each GTEx tissue. These databases have been made available for download at PredictDB.org.

99Supplementary Data 2 shows the list of GWA/GWAMA studies we considered in this analysis. We applied S-PrediXcan to these studies using the transcriptome models trained on GTEx studies for patched version 6. For simplicity, S-PrediXcan only considers those SNPs that have a matching set of alleles in the prediction model, and adjusts the dosages (2 dosage) if the alleles are swapped.

100To make the results of this study broadly accessible, we built a Postgre SQL relational database to store S-PrediXcan results, and serve them via a web application http://gene2pheno.org.

101We also applied SMR to the same set of GWAMA studies, using the GTEx eQTL associations. We downloaded version 0.66 of the software from the SMR website, and ran it using the default parameters. We converted the GWAMA and GTEx eQTL studies to SMR input formats. In order to have SMR compute the colocalization test, for those few GWAMA studies where allele frequency was not reported, we filled in with frequencies from the 1000 Genomes Project as an approximation. We also used the 1000 Genomes genotype data as reference panel for SMR.

102Next we ran COLOC (as downloaded from the Comprehensive R Archive Network) over the same set of GWAMA and eQTL studies. We used the Approximate Bayes Factor colocalization analysis, with effect sizes, their standard errors, allele frequencies and sample sizes as arguments. When the frequency information was missing from the GWAS, we filled in with data from the 1000 Genomes Project.

103For comparison purposes, we have also included the results of the application of Summary-TWAS to 30 traits publicly shared by the authors.

104Formal similarity with TWAS can be made more explicit by rewriting S-PrediXcan formula in matrix form. With the following notation and definitionsand correlation matrix of SNPs in the model for gene g

105it is quite straightforward to write the numerator in (1) and (11) as

106and in the denominator, the variance of the predicted expression level of gene g, as

107Thus

108This equation has the same form as the TWAS expression if we use the scaled weight vector instead of Wg. Summary-TWAS imputes the Z-score for the gene-level result assuming that under the null hypothesis, the Z-scores are normally distributed with the same correlation structure as the SNPs; whereas in S-PrediXcan we compute the results of PrediXcan using summary statistics. Thus, S-TWAS and S-PrediXcan yield equivalent mathematical expressions (after setting the factor ).

109The S-PrediXcan formula when only the top eQTL is used to predict the expression level of a gene can be expressed aswhere Z1 is the GWAS Z-score of the top eQTL in the model for gene. Thus

110SMR quantifies the strength of the association between expression levels of a gene and complex traits with TSMR using the following function of the eQTL and GWAS Z-score statistics

111Here is the Z-score (=effect size/standard error) of the association between SNP and gene expression, and is the Z-score of the association between SNP and trait.

112This SMR statistic (TSMR) is not a random variable as assumed in ref.. To prove this, we performed simulations following those described in ref.. We generated 105 pairs of values for and. was sampled from a distribution. was sampled from a non-central distribution with parameter λ = 29 (a value chosen to mimic results from, see ref. ). Only pairs with eQTLs satisfying genome-wide significance (p < 5 × 10−8) were kept. We performed a QQ plot and observed deflation when comparing to random values sampled from a distribution (Fig. 5e). This simulation was repeated 1000 times, and we observed a mean of TSMR close to 0.93 (Fig. 5f).

113Only in two extreme cases, the chi-square approximation holds, when or. In these extremes, we can apply Taylor expansions to find an interpretable form of the SMR statistic.

114If, i.e., the eQTL association is much more significant than the GWAS association,

115so that for large enough relative to,

116using Eq. 12. Thus, in this case, the SMR statistic is slightly smaller than the (top eQTL based) S-PrediXcan -square. This reduced significance is accounting for the uncertainty in the eQTL association. As the evidence for eQTL association grows, the denominator increases and the difference tends to 0.

117On the other extreme when the GWAS association is much stronger than the eQTLs,,

118so that analogously

119In both extremes, the SMR statistic significance is approximately equal to the less significant of the two statistics GWAS or eQTL, albeit strictly smaller.

120In between the two extremes, the right distribution must be computed using numerical methods. When we look at the empirical distribution of the SMR statistic’s p-value against the GWAS and eQTL (top eQTL for the gene) p-values, we find the ceiling of the SMR statistic is maintained as shown in Fig. 5e, f.

121Genotype files were obtained from dbGaP, and updated to release 35 of the probe annotations published by Affymetrix via PLINK. Probes were filtered out that had a minor allele frequency of <0.01, were missing in >10% of subjects, or did not fit Hardy-Weinberg equilibrium. Subjects were dropped that had an unexpected level of heterozygosity (F>0.05). Finally the HRC-1000G-check-bim.pl script (http://www.well.ox.ac.uk/~wrayner/tools/) was used to perform some final filtering and split data by chromosome. Phasing (via eagle v2.3) and imputation against the HRC r1.1 2016 panel (via minimac3) were carried out by the Michigan Imputation Server.

122European samples had been split into ten groups during imputation to ease the computational burden on the Michigan server, so after obtaining the imputed.vcf files, we used the software PLINK to convert the genotype files into the PLINK binary file format and merge the ten groups of samples together, while dropping any variants not found in all sample groups. For the association analysis, we performed a logistic regression using PLINK, and following QC practices from ref. we filtered out individuals with genotype missingness >0.03 and filtered out variants with minor allele frequency <0.01, missingness >0.05, out of Hardy-Weinberg equilibrium significant at 1e-6, or had imputation quality <0.8. We used gender and the first ten genetic principal components as obtained from dbGaP as covariates. Following all filtering, our analysis included 61,444 European samples with 7,120,064 variants. MetaXcan was then applied to these GWAS results, using the 45 prediction models (GTEx and DGN).

123We make our software publicly available on a GitHub repository: https://github.com/hakyimlab/MetaXcan. A short working example can be found on the GitHub page; more extensive documentation can be found on the project’s wiki page.

124The underlying GWAS results used in this analysis were downloaded from publicly available resources listed in Supplementary Data 2. The relevant GTEx gene expression data was obtained from dbGAP using accession phs000424.v6.p1. The GERA study was downloaded from dbGAP using accession number phs000674.v2.p2. WTCCC data was downloaded from WTCCC EGA european genome-phenome archive.

125The list of ClinVar genes was downloaded from https://www.ncbi.nlm.nih.gov/clinvar/. TWAS results published in ref. were used. Prediction model weights and covariances for different tissues are available from the predictdb.org resource. The results of MetaXcan applied to the 44 human tissues and a broad set of phenotypes can be queried in gene2pheno.org, and we make the full data set of results available via the public GitHub repository https://github.com/hakyimlab/MetaXcan.