PrediXcan: Trait Mapping Using Human Transcriptome Regulation

Genome-wide association studies (GWAS) have identified thousands of variants robustly associated with complex traits. However, the biological mechanisms underlying these associations are, in general, not well understood. We propose a gene-based association method called PrediXcan that directly tests the molecular mechanisms through which genetic variation affects phenotype. The approach estimates the component of gene expression determined by an individual’s genetic profile and correlates the “imputed” gene expression with the phenotype under investigation to identify genes involved in the etiology of the phenotype. The genetically regulated gene expression is estimated using whole-genome tissue-dependent prediction models trained with reference transcriptome datasets. PrediXcan enjoys the benefits of gene-based approaches such as reduced multiple testing burden, more comprehensive annotation of gene function compared to that derived from single variants, and a principled approach to the design of follow-up experiments while also integrating knowledge of regulatory function. Since no actual expression data are used in the analysis of GWAS data - only in silico expression - reverse causality problems are largely avoided. PrediXcan harnesses reference transcriptome data for disease mapping studies. Our results demonstrate that PrediXcan can detect known and novel genes associated with disease traits and provide insights into the mechanism of these associations.


INTRODUCTION
Genome-wide association studies (GWAS) have been remarkably successful in identifying susceptibility loci for complex diseases. These studies typically conduct single-variant tests of association to interrogate the genome in an agnostic fashion and, due to modest effect sizes, have come to rely on ever-greater sample sizes 1,2 to make meaningful inferences. We have been less successful in developing methods that improve on existing simple approaches. In general, the genetic associations identified as genome-wide significant thus far account for only a modest proportion of variance in disease risk 3 . Indeed, there is now widespread recognition, if not consensus, that GWAS of disease susceptibility (for which, the relevant genetic effects may be very small) and pharmacologic traits (for which large effect sizes are not unusual) 4,5 have resulted in limited conclusive findings on the genetic factors contributing to complex traits. Importantly, the functional significance of most discovered loci, including even those that have been the most reproducibly associated, remains unclear. Assigning a causal link to the nearest gene falls short of elucidating a functional connection, as recently demonstrated by the obesity-associated variants within FTO that form longrange functional connections with IRX3 6 . And while GWAS will no doubt continue to identify many more susceptibility loci, the question of how to advance biological knowledge of the underlying mechanisms of disease risk remains a paramount challenge.
A large portion of phenotypic variability in disease risk for a broad spectrum of disease phenotypes can be explained by regulatory variants, i.e. genetic variants that regulate the expression levels of genes 7-10 . For example, almost 80% of the chip-based heritability of disease risk for 11 diseases from the WTCCC can be explained by genome variation in DNase I hypersensitivity sites, which are likely to regulate chromatin accessibility and thus transcription 11 . Large genomic consortia (e.g., ENCODE 12 ) are generating an unprecedented volume of data on the function of genetic variation. The Genotype-Tissue Expression (GTEx 13 ) project is an NIH Common Fund project that aims to collect a comprehensive set of tissues from 900 deceased donors (for a total of about 20,000 samples) and to Methodologically, gene-based approaches and multi-marker association tests have been developed as alternatives to traditional single-variant tests. By conducting tests of association on biologically informed aggregates of SNPs, such tests seek to evaluate a priori functionally relevant units of the genome and, in many cases, reduce the multiple-testing penalty that plague single-variant approaches, by 10 to 100 fold.
The incorporation of -omics data, such as those being generated by high-resolution transcriptome studies, provides a means to extend genome-wide association studies by addressing the functional gap. Technological advances in high-throughput methods have reinforced the important finding that intermediate molecular phenotypes are under significant genetic regulation, with expression quantitative trait loci (eQTLs) as the predominant example. However, approaches that fully leverage the comprehensive regulatory knowledge generated by transcriptome studies are relatively lacking despite the fact that these studies have the potential to dramatically improve our understanding of the genetic basis of complex traits 13 .
We hypothesized that a SNP aggregation approach that integrates information on whether a SNP regulates the expression of a gene can greatly increase the power to identify trait-associated loci either from a strong functional SNP signal or from a combination of modest signals, the so-called grey area of GWAS. The present study suggests that PrediXcan, a novel method that incorporates information on gene regulation from a set of markers, increases the power to detect associations relative to traditional SNP-based GWAS and known gene-based tests under a broad range of genetic architectures and provides mechanistic insights and more easily interpreted direction of effect into the observed associations.

PrediXcan method
PrediXcan, by design, exploits genetic control of phenotype through the mechanism of gene regulation as a way to identify trait-associated genes. Figure 1 is a schematic diagram of the regulatory mechanism that is tested with PrediXcan. An individual's gene expression level (typically unobserved in a GWAS) is decomposed into a genetically regulated expression (GReX) component, a component altered by the trait itself (i.e., a reverse causal effect that may occur if disease status or other conditions alter expression levels), and the remaining component attributable to environmental and other factors. PrediXcan tests the mediating effect of gene expression by quantifying the association between GReX and the phenotype of interest.
We use reference transcriptome datasets from studies such as GTEx 13 , GEUVADIS 15 , and DGN 16 among others, to train additive models of gene expression levels. These models allow us to estimate the genetically regulated expression, GReX.
We denote the estimated value with a hat, . These estimates constitute multiple-SNP prediction of expression levels. The weights for the estimation are stored in our publicly available database.
The analogy with genotype imputation is relevant here. Genotype imputation uses information from a reference sample to learn how to impute genotypes at the unmeasured SNPs in the test set. Similarly, PrediXcan uses a reference dataset in which both genome variation and gene expression levels have been measured to develop prediction models for gene expression. We use these prediction models to "impute" gene expression (which is unobserved in a typical GWAS), and we do so by estimating the genetically determined component, GReX.
PrediXcan application to a GWAS dataset consists of "imputing" the transcriptome using the weights derived from reference transcriptome datasets and correlating the GReX with the phenotype of interest using regression methods (e.g., linear, logistic, Cox) or non-parametric approaches (e.g., Spearman). (For the specific results on disease phenotypes analyzed here, we used logistic regression with disease status.) We are aware of the attenuation bias that arises because of the error in the estimation of GReX. This is a subject to be investigated in the future, but this bias does not invalidate our analysis since we only use the estimate of GReX as a discovery tool. Figure 2 summarizes the flow of the method development described above.

Features of PrediXcan
PrediXcan is, as we have emphasized, particularly focused on a mechanism -gene expression regulation -that has already been established as being contributory to common diseases, including psychiatric and neurodevelopmental disorders 7 . The test has the potential to identify gene targets for therapeutic applications because it is inherently mechanism-based and provides directionality.
Additional advantages include: • Like other gene-based tests, it has much smaller multiple-testing burden (~20K tests maximum, ~10K genes with high quality prediction in most tissues) compared to single variant tests (~5-10M tests). Moving beyond the stringent Bonferroni correction, priors on genes can be less restrictive than for SNPs.
• Informative priors and groupings of functional units (based on known pathways, for example) are much more straightforward to construct for genes than SNPs.
• No actual transcriptome data are required since the predicted expression levels are a function of genetic variation alone. Thus, the method can be applied to any existing dataset with large-scale genome interrogation such as those in dbGaP or other repositories. Re-analyses of existing datasets, with a focus on mechanism using PrediXcan, address a gap that has largely characterized GWAS to date.
• Reverse causality is not a major concern since disease status or drug treatment does not alter germline genomic variation.
• Meta-analysis of gene-based results is simplified since less stringent harmonization between studies is required.
• Multiple tissues can be evaluated using a reference transcriptome dataset (such as GTEx). In general, the only limitation is the availability of gene expression data in the given tissue for model building, which need not be, from the same study as that used for phenotype investigation. In cases where transcriptome data are available, separate analyses should be performed to simplify interpretation.
• The approach can be applied to common or rare variants. In general, larger sample sizes for the training set will be needed to achieve good prediction models with rare variants.

Publicly available database of prediction models and software
We make the prediction models (derived from LASSO 18 and elastic net 19 ) and the software to predict the transcriptome (in a variety of tissues) (see Materials and Methods) publicly available (see https://github.com/hakyimlab/PrediXcan).

Prediction model selection
We built prediction models in the DGN whole blood cohort using LASSO, the elastic net (⍺=0.5), and polygenic score at several p-value thresholds (single top SNP, 1e-04, 0.001, 0.01, 0.05, 0.5, 1). We assessed predictive performance using 10-fold cross-validation (R 2 of estimated GReX vs. observed expression) as well as in an independent set. We found that LASSO performed similarly to the elastic net and that LASSO outperformed the polygenic score at all thresholds, although all methods are highly correlated (see Supplemental Figure 1). For subsequent analyses, we focused on the prediction models using the elastic net because we found it to perform well and to be more robust to slight changes in input SNPs (potentially due to variations in imputation quality between cohorts).

Heritability estimation and comparison with prediction R 2
We estimated the heritability of gene expression in DGN attributable to SNPs in the vicinity of each gene using a mixed-effects model (see Materials and Methods) and calculated variances using restricted maximum likelihood as implemented in GCTA 20 .
We use only local SNPs since we found that heritability estimates using all genotyped SNPs were too noisy to make meaningful inferences.
We use heritability estimates as our benchmark for the prediction R 2 since this constitutes the upper limit of our prediction performance. For genes for which an elastic net model was available (n=10,427), the average heritability in DGN was 0.153. In comparison, the average 10-fold cross-validated prediction R 2 for elastic net was quite close at 0.137; for the polygenic score (P<1x10 -4 ) and top-SNP models, average prediction R 2 values were sizably lower at 0.099 and 0.114, respectively. We show the performance R 2 for each model in Figure 3, with the corresponding heritability estimate and confidence interval in the background for comparison. We also note that elastic net predictive performance reached or exceeded the lower bound of the heritability estimate for 94% of genes, while polygenic score (P<1x10 -4 ) did so for just 76% of the genes and the top SNP for 80% of the genes (Figure 3), consistent with the performance ranking given by the average (across genes) R 2 .

Comparing imputed vs. directly genotyped SNPs as reference set for prediction model
Predictive performance of elastic net was similar whether all SNPs from the 1000 Genomes imputation or the HapMap Phase II subset were included in the model building (Supplemental Figure 2). Models based on imputed data (both the 1000 Genomes and the HapMap subset) substantially outperformed models based on genotyped SNPs in WTCCC (Supplemental Figure 2). Thus, we chose the elastic net models built in the smaller HapMap SNP subset, relative to 1000 Genomes, in our applications of PrediXcan to reduce computation time without sacrificing performance.
As reference transcriptome studies increase in sample size, we may need to switch to a more dense imputation to take advantage of increased prediction performance from rare variants.

Prediction performance in an independent test set
We also tested the prediction models trained in the DGN whole blood cohort on several independent test cohorts with available whole-genome genotype and transcriptome data. We used weights derived from the DGN whole blood data ("training set") to predict gene expression levels (treated as quantitative traits) in GEUVADIS LCLs (lymphoblastoid cell lines) and nine GTEx pilot tissues ("test sets"). Figure 4 provides a Q-Q plot showing the expected (under the null, correlation between two independent vectors with the same sample size) and observed R 2 (between observed and predicted) from the elastic net prediction in GEUVADIS LCLs. We find a substantial (muscle), 0.0422 (tibial nerve), 0.0374 (sun exposed skin), 0.0398 (thyroid), and 0.0458 (whole blood). Interestingly, we also find a substantial departure from the null distribution of expected R 2 values for predicted expression using DGN weights in each of the nine GTEx tissues suggesting that models developed in whole blood are still useful for understanding diseases that affect other primary tissues (Supplemental Figure 3). Consistent with this, average prediction R 2 is highest for whole blood as expected but the loss in power for other tissues is modest. Figure 5 illustrates the genes with some of the highest correlations from this analysis, providing a comparison of the predicted expression and the observed expression. Among these genes, both ERAP2 and its paralog ERAP1 play fundamental roles in MHC antigen presentation 21 , immune activation and inflammation.

Marginal gain in performance when adding distant SNPs
We also generated prediction models trained in the DGN whole blood cohort that included trans-eQTLs (>1Mb from gene start or end or on a different chromosome) generated from linear regression (p<10 -5 ). We tested the predictive performance of these models in the GTEx whole blood cohort. While a few genes had higher correlations between predicted and observed expression than expected by chance, the departure from the null distribution was much smaller than that for the prediction models based on local SNPs (Supplemental Figure 4)

Comparison of gene-based tests
We applied PrediXcan and two widely-used gene-based tests (VEGAS and SKAT) to WTCCC. In a Q-Q plot showing all three distributions of p-values, for genes outside of the HLA region, from these gene-based tests (Figure 7), SKAT had improved performance relative to VEGAS, and PrediXcan showed the most extreme departure from the null at the tail end of the distribution.

Replication of PrediXcan findings
To replicate our findings, we applied the DGN elastic net whole blood prediction models to an independent rheumatoid arthritis GWAS from Vanderbilt University's BioVU repository (see Materials and Methods). Both genes (DCLRE1B and PTPN22) that were found to be genome-wide significant in the WTCCC rheumatoid arthritis data were also significant, with concordant direction of effect, in the replication samples (p = 0.012 and p = 0.036, respectively). Application of the method to WTCCC data recapitulated many known loci but also identified novel genome-wide significant genes. We believe that a systematic reanalyses of GWAS datasets in comprehensive repositories such as dbGAP and the European Genome-phenome Archive (EGA) could provide a cost-effective approach to uncovering novel disease mechanisms using only existing genomic resources.

DISCUSSION
In contrast to other gene-based tests, PrediXcan provides the direction of effect, which may yield opportunities for therapeutic development. The development of therapeutics that down-regulate a gene is generally easier to achieve than therapeutics that up-regulate a gene; thus, genes with expression levels that are positively correlated with disease risk may be more favorable drug targets for novel therapies. The direction of effect may also provide information to elucidate pathways and the opportunity to explore systems-based approaches to the development of disease. The prediction models can be applied to genotype data of subjects in large biobanks to investigate potential side-effects of drugs with specific gene targets. Finally, direction of effect can be used to improve the interpretation of sequence analyses of genes showing significant correlation of predicted expression with phenotype, since phenotypes associated with reduced expression of genes are more likely to show a relative excess of rare variants. Indeed, we believe that PrediXcan offers intriguing opportunities to combine results of rare and common variant association tests within whole genome sequencing studies, and more generally, to combine results of rare variant gene-based tests from sequencing studies with results of PrediXcan gene-based tests from the large body of existing GWAS for the same phenotypes. Thus, PrediXcan is a method developed to integrate -omics data that can facilitate integration of results from common and rare variant studies.
Regarding the multiple testing correction approach, here we have used Bonferroni correction using the total number of genes tested. In general, both singlevariant and PrediXcan analyses will be performed; thus the question that arises is how to address the issue of multiple testing adjustment. The prior probability for a SNP to be causal is much smaller than the prior probability of causality for a gene so it would not be fair to subject SNP tests and gene-based tests to the same level of adjustment.
Since we are presenting only gene-based results in our application and given the highly conservative nature of Bonferroni correction, there is no need to further adjust our results. A more conservative approach would be to divide the significance threshold used by a factor of two for the multiple testing using gene-based and SNP-based approaches.
Given the large contribution of regulatory variants on complex traits 9,10,37 , our method is likely to identify causal genes. However, we do not claim causality since SNPs that contribute to the expression of a gene can also act through other mechanisms to determine the phenotype of interest. Replication and experimental validations are needed to determine causality.
In conclusion, we presented a novel gene-based test, PrediXcan that incorporates functional information with regard to gene regulation to identify genes associated with disease traits in large GWAS or whole genome sequence datasets. Our method has the advantage of providing biological insights into the mechanism, namely regulation of gene expression, and direction of effect. This approach can be readily applied to existing GWAS datasets through the use of our publically available PredictDB resource. We further show the utility of our approach by identifying and replicating a number of novel candidate associations within the previously analyzed WTCCC dataset.

DGN RNA-Seq Dataset
We obtained whole blood RNA-Seq 38 and genome-wide genotype data for 922 individuals from the Depression Genes and Networks cohort 16 , all of European ancestry.
For our analyses, we used the HCP (hidden covariates with prior) normalized gene-level expression data used for the trans-eQTL analysis in Battle et al. 16

GEUVADIS RNA-Seq Dataset
We obtained freely available RNA-Seq data from 421 lymphoblastoid cell lines (LCLs) generated by the GEUVADIS consortium 15 and genotype data generated by the 1000 Genomes project (http://www.geuvadis.org/web/geuvadis/RNAseq-project). We used GEUVADIS as a validation dataset to test the gene prediction models generated in the DGN cohort.

GTEx RNA-Seq Datasets
We used the nine tissues with the largest sample size in the Genotype-Tissue Expression (GTEx) Pilot Project 14 to test the gene prediction models generated in the DGN cohort. Tissue samples included subcutaneous adipose (n=115), tibial artery (n=122), left ventricle heart (n=88), lung (n=126), skeletal muscle (n=143), tibial nerve (n=98), skin from the sun-exposed portion of the lower leg (n=114), thyroid (n=112), and whole blood (n=162). In each tissue, normalized gene expression was adjusted for gender, the top 3 principal components (derived from genotype data), and the top 15 PEER factors (to quantify batch effects and experimental confounders) 41 . We used GTEx to test the portability of predictors developed in whole blood (from the DGN cohort) across a wide variety of tissues.

Additive model for gene expression traits
We use an additive genetic model to characterize gene expression traits: where Y g is the expression trait of gene g, w k,g is the effect size of marker k for gene g, X k is the number of reference alleles of marker k, and ϵ is the contribution of other factors that determine the expression trait assumed to be independent of the genetic component. We note that the summation in model (1)  in DGN shared genetic relatedness ( ) in excess of 5% and thus all were included in the narrow-sense heritability estimation. SNPs in the vicinity of each gene (within 1Mb of gene start or end, as defined by the GENCODE 45 version 12 gene annotation), with MAF > 0.05, and in Hardy-Weinberg Equilibrium (P > 0.05) were used to construct the GRM for each gene. We calculated the proportion of the variance of gene expression explained by these local SNPs using the following mixed-effects model 37 : where Y is a gene expression trait and b a vector of fixed effects. Here is the GRM calculated from the local SNPs, and (the random effect) denotes the genetic effect attributable to the set of local SNPs with var ( ) = .
In this paper we focus on the component of heritability driven by SNPs in the vicinity of each gene since the component based on distal SNPs could not be estimated with enough accuracy to make meaningful inferences.

Estimation of the genetic component of gene expression levels (GReX)
In the simple polygenic score approach, we estimate w k as the single-variant coefficient derived from regressing the gene expression trait Y on variant X k (as implemented in the eQTL analysis software Matrix eQTL 46 ) using the reference transcriptome data. This yields an estimate, , for a GWAS study sample, of the (unobserved) genetically determined expression of each gene g: In this implementation of polygenic score, we include all SNPs (regardless of linkage disequilibrium [LD]) that are associated with the expression level of the gene at a chosen p-value threshold in the prediction model.
In contrast, LASSO uses an L1 penalty as a variable selection method to select a sparse set of (uncorrelated) predictors 18 while the elastic net linearly combines the L1 and L2 penalties of LASSO and ridge regression respectively to perform variable selection 19 . We used the R package glmnet to implement LASSO and elastic net with ⍺=0.5 (http://www.jstatsoft.org/v33/i01).
For each gene, LASSO, the elastic net and the simple polygenic score were used to provide an estimate of (using equation 2, with the effect size estimates , , , and , , respectively). We included only local SNPs (within 1Mb of the gene start or end). In order to determine the optimal modeling method, we compared the 10-fold cross-validated prediction R 2 (the square of the correlation between predicted and observed expression) for the simple polygenic score ( ) at several p-value thresholds (single top SNP, 1x10 -4 , 0.001, 0.01, 0.05, 0.5, 1) with that from LASSO ( ) and elastic net ( ).
We also compared the 10-fold cross-validated prediction R 2 from elastic net models with different starting SNP sets from the DGN genotype imputation (4.6M 1000 Genomes Project SNPs (MAF>0.05, R 2 >0.8, non-ambiguous strand), the 1.9M of these SNPs that are also in HapMap Phase II, and the 300K of these SNPs that were genotyped in the WTCCC).

Performance of transcriptome prediction in independent cohorts
We tested the feasibility of predicting the transcriptome (i.e., estimating the genetic component of each gene expression trait, , in an independent test transcriptome dataset) using the elastic net effect sizes trained in the DGN whole blood data (n=922). For the test sets, we used independent RNA-Seq datasets from 421 LCL cell lines from the 1000 Genomes project generated by the GEUVADIS consortium 15 and the nine tissues from the GTEx pilot project 14 (see Supplemental Figure 3). To assess performance, we used the square of the Pearson correlation, R 2 , between predicted and observed expression levels.

PrediXcan in the WTCCC GWAS Datasets
To illustrate the method, we applied gene prediction models (derived from whole blood) consisting of DGN elastic net predictors to the seven WTCCC disease studies --bipolar disorder (BD), coronary artery disease (CAD), hypertension (HT), type 1 diabetes (T1D), type 2 diabetes (T2D), Crohn's disease (CD), and rheumatoid arthritis (RA) 22 . Genotypes imputed to the 1000G reference sets were used. Imputation was done using the University of Michigan Imputation-Server and the same parameters as described for the imputation of DGN data. For each disease, cases and controls (1958 Birth Cohort and the UK Blood Service Cohort) were jointly imputed to avoid subtle differences between cases and controls not attributable to disease risk. We excluded all SNPs with an imputation R 2 < 0.8 and for computational speed we kept only the HapMap Phase II subset of SNPs.
For each WTCCC disease, we estimated , and tested it for association with disease risk using logistic regression in R (R-project.org). We restricted our PrediXcan analysis to include genes with a cross-validated prediction R 2 > 0.01 (10% correlation) in the DGN sample. Because the WTCCC studies use shared controls, pleiotropy analyses using these datasets would not be straightforward, and comparison of results across diseases was avoided.

GWAS Enrichment analysis
Relative to recent association studies, the WTCCC has a small sample size (~2,000 cases and ~3,000 controls per disease). Thus, even with our novel method and a reduced multiple testing burden, our ability to detect numerous novel gene associations may be limited. Alternatively, we tested each disease for an enrichment of known disease genes identified from the NHGRI GWAS Catalog 25 . For each disease, we used the reported genes from the GWAS catalog as the set of known disease genes. We excluded studies listed in the NHGRI GWAS catalog that included the WTCCC samples in order to make sure our known gene lists were independent from the current analysis. We then counted the number of known disease genes that had a PrediXcan p-value below a given threshold. We compared this to the null expectation based on 10,000 randomly drawn gene sets of similar size to the known disease gene set to derive an enrichment p-value. We tested enrichment using PrediXcan p-value thresholds of 0.05 and 0.01.

Comparison to large single variant meta-analyses
For the top PrediXcan results in the WTCCC, we cross-referenced the SNPs in the prediction models for these genes with the publically available single-SNP metaanalysis summary results. We excluded T1D from this analysis because, to our knowledge, there are no publically available meta-analysis studies of this disease. We used meta-analyses results for systolic and diastolic blood pressure as a proxy for hypertension. For CD, RA, and BD we were able to use meta-analyses for the same diseases (CD 47 : http://www.ibdgenetics.org/downloads.html, RA 48 , and BD 31 : http://www.med.unc.edu/pgc/downloads).

Comparison of gene-based tests (PrediXcan, SKAT, VEGAS)
We compared the results derived from PrediXcan with those from two widelyused gene-based tests, namely VEGAS 49 and SKAT 50,51 . VEGAS aggregates information from the full set of SNPs within a gene and accounts for LD using simulations from the multivariate normal distribution. SKAT is a kernel-based association test that evaluates the regression coefficients of the SNPs within a gene by a variance component score test in a mixed model framework. We generated BEDformatted files for SNPs and genes (as defined by GENCODE v12) and mapped SNPs that met post-imputation QC parameters to gene regions using bedtools. The use of an offline Perl implementation for VEGAS allowed us to examine the dependence of the results from this approach on LD information through the use of the actual genotype data (versus the default HapMap CEU reference panel data). We developed an Rbased pipeline that invokes the SKAT package (version 1.0.1) that is publicly available from CRAN. We generated a Q-Q plot showing the distribution of gene-level p-values for association with RA (for genes outside the HLA region) derived from each genebased test to test for systematic departure from the null expectation (of uniform pvalues).

Replication of PrediXcan findings
We   The sample size of the study is denoted by n, m is the number of genes considered, M is the total number of SNPs, and p is the number of available tissues. The second panel shows the additive model used to build a database of prediction models, PredictDB. T represents the expression trait, and X k is the number of reference alleles for SNP k. The coefficients of the models for each tissue are fitted using the reference transcriptome datasets and optimal statistical learning methods chosen among LASSO, Elastic Net, OmicKriging, etc. The bottom panel shows the application of PrediXcan to a GWAS dataset. Using genetic variation data from the GWAS and weights in PredictDB, we "impute" expression levels for the whole transcriptome. These imputed levels are correlated with the trait using regression (e.g., linear, logistic, Cox) or non-parametric (Spearman) approaches. (For the disease phenotypes in the WTCCC datasets and the replication dataset reported here, we used logistic regression with disease status.) Figure 3. Cross-validated prediction performance vs heritability. This figure shows the prediction performance (R2 of GReX vs. observed expression in red) compared to gene expression heritability estimates (black with 95% confidence interval in gray). Performance was assessed using 10-fold cross-validation in the DGN whole blood cohort (n=922) with the elastic net, polygenic score (p < 1e-04), and using the top SNP for prediction.

Figure 4. Prediction performance of elastic net tested on a separate cohort.
Using whole blood prediction models trained in DGN, we compared predicted levels of expression with observed levels on lymphoblastoid cell lines from the 1000 Genomes project. RNA-sequenced data (n=421) on these cell lines have been made publicly available by the GEUVADIS consortium. Left panel shows the squared correlation, 2 , between predicted and observed levels plotted against the null distribution of 2 . Right panel shows prediction performance (R 2 of GReX vs. observed expression in green) compared to GEUVADIS gene expression heritability (h 2 ) estimates (black with 95% confidence interval in gray).

Figure 5. Examples of well-predicted genes.
These plots show observed vs. predicted levels of 4 genes. Predicted levels were computed using whole blood elastic net prediction models trained in DGN data. Observed levels were RNA-seq data in lymphoblastoid cell lines generated by the GEUVADIS consortium.   Chromosome and gene start position are based on GENCODE version 12. The cross validated prediction R 2 between predicted and observed gene expression is based on 10-fold cross validation within the DGN whole blood sample. SNPs listed below are included in the DGN whole blood prediction model for PTPRE.
The association p-values for these SNPs in the WTCCC data and the PGC Bipolar disorder meta-analysis are shown.
Supplemental Figure 1. Comparison of 10-fold cross-validated predictive performance between all tested methods (LASSO, elastic net with α=0.5, top SNP, polygenic score at several p-value thresholds) in the DGN whole blood cohort. Predictive performance was measured by the R 2 between predicted (GReX) and observed expression.
Supplemental Figure 2. Comparison of 10-fold cross-validated predictive performance of elastic net in different starting SNP sets (4.6M 1000 Genomes Project (TGP) SNPs, 1.9 M HapMap Phase II SNPs, 300K WTCCC genotyped SNPs) in the DGN whole blood cohort. Predictive performance was measured by the R 2 between predicted (GReX) and observed expression.

Supplemental Figure 3. Prediction performance of elastic net in GTEx tissues.
Using whole blood prediction models trained in DGN, we compared predicted levels of expression with observed levels from nine tissues of the GTEx pilot project. The observed squared correlation between predicted and observed gene expression levels, 2 , is plotted against the null distribution of 2 .

Supplemental Figure 4. Comparison of prediction performance between localand distal-based prediction models.
Using whole blood prediction models trained in DGN, we compared predicted levels of expression with observed levels in GTEx whole blood. Local predictors were generated using elastic net on SNPs within 1Mb of each gene and distal predictors include any trans-eQTLs outside this region with a linear regression p<10 -5 . The observed (y-axis) squared correlation between predicted and observed gene expression levels, 2 , is plotted against the null distribution of 2 (x-axis).