Identifying gene targets for brain-related traits using transcriptomic and methylomic data from blood

ARTICLE Identifying gene targets for brain-related traits using transcriptomic and methylomic data from blood Ting Qi et al. # Understanding the difference in genetic regulation of gene expression between brain and blood is important for discovering genes for brain-related traits and disorders. Here, we estimate the correlation of genetic effects at the top-associated cis-expression or -DNA methylation (DNAm) quantitative trait loci (cis-eQTLs or cis-mQTLs) between brain and blood (r b ). Using publicly available data, we find that genetic effects at the top cis-eQTLs or mQTLs are highly correlated between independent brain and blood samples (^ rb¼0:70 for ciseQTLs and ^ rb¼0:78 for cis-mQTLs). Using meta-analyzed brain cis-eQTL/mQTL data (n= 526 to 1194), we identify 61 genes and 167 DNAm sites associated with four brain-related phenotypes, most of which are a subset of the discoveries (97 genes and 295 DNAm sites) using data from blood with larger sample sizes (n=1980 to 14,115). Our results demonstrate the gain of power in gene discovery for brain-related phenotypes using blood cis-eQTL/mQTL data with large sample sizes. DOI: 10.1038/s41467-018-04558-1 OPEN Correspondence and requests for materials should be addressed to J.Y. (email: [email protected]). #A full list of consortia members appears at the end of the paper. NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications 1 1234567890():,; Genome-wide association studies (GWAS) have discovered thousands of genetic variants associated with complex traits and diseases1–3. Most trait-associated variants reside in non-coding regions of the genome4,5, suggesting that genetic variants may affect the trait through regulation of gene expression6,7. With the advances in microarray and sequencing technologies, genome-wide genotype and gene expression data available from relatively large samples have been generated to identify genetic variants affecting transcription abundance8–10, i.e., expression quantitative trait loci (eQTLs). Current eQTL studies are biased toward the most accessible tissues (e.g., blood), which are often not the most relevant tissues to the traits and diseases of interest. The Genotype-Tissue Expression (GTEx) project11–13 provides a comprehensive resource of data to investigate the genetic variation of gene expression across a broad range of tissues and cell types. Recent studies have utilized the GTEx data to demonstrate that genetic correlation of gene expression between tissues in local regions (i.e., ±1Mb of the transcription start site) is much higher than that in distal regions14, consistent with the conclusions from the latest GTEx release13, and that there is no evidence for the tissue-relevant eQTLs being enriched for associations with complex traits15. For studies that integrate GWAS results with eQTL or DNA methylation QTL (mQTL) data to identify putative functional genes and regulatory elements for brain-related phenotypes and diseases16,17, the statistical power is limited by the small sample sizes of the brain eQTL or mQTL data (often in the order of 100s). On the other hand, there are blood eQTL and mQTL data available from thousands of individuals8,9and the sample sizes of some of the ongoing projects have reached 10,000s (e.g., the GoDMC and eQTLGen consortia). The questions are to what extent the cis-genetic effects on gene expression and DNA methylation (DNAm) in blood differ from those in brain and whether we can gain power to detect associations of genes (or DNAm sites) with brain-related traits by using the cis-eQTL (or cis-mQTL) effects estimated from a large blood sample as proxies for those in brain. Liu et al.14 extended the stratified linkage disequilibrium (LD) score regression method to estimate genetic correlation (r g ) of gene expression between tissues at all SNPs in local or distal regions and showed that the mean estimate of pairwise r g at all local SNPs (i.e. cis-genetic correlation) was ~0.75 in 11 GTEx tissues but they did not estimate r g between brain and blood. In this study, we use a summary-data-based method to estimate the correlation of effect sizes of the top-associated ciseQTLs (or cis-mQTLs) between blood and brain for genes expressed (or CpG sites methylated) in both tissues, accounting for errors in their estimated effects. We demonstrate by simulation and analysis of real data the gain of power by using cis-eQTL or cis-mQTL effects estimated in blood as proxies of those in brain to identify putative functional genes for brain-related complex traits and diseases. Results Correlation of cis-eQTL effects between brain and blood.To quantify the similarity of genetic effects at the top-associated ciseQTLs (or cis-mQTLs) between two tissues, we used a summarydata-based approach to estimate the correlation of cis-effects between two tissues (r b ) correcting for errors in the estimated ciseQTL (or cis-mQTL) effects and sample overlap (Supplementary Fig. 1and Methods). We showed by simulation (Supplementary Note 1) that r b is a good estimator of correlation of the true values of cis-genetic effects (Supplementary Fig. 2). Note that the r b method is distinct from the Spearman or Pearson correlation approach13 because the latter does not account for errors in the estimated eQTL effects and thereby leads to an underestimation of the correlation of true eQTL effects. We applied our method to estimate ^ rbat the top cis-eQTLs between different brain regions and between brain and blood in one data set, and between brain and blood in two data sets using summary-level data from GTEx v6 (whole blood and 10 brain regions)11, the CommonMind Consortium (CMC; dorsolateral prefrontal cortex)18, the Religious Orders Study and Memory and Aging Project (ROSMAP) 19, and the Brain eQTL Almanac project (Braineac; 10 brain regions)20 (Methods and Supplementary Table 1). All eQTL effects were re-scaled based on the expression level per gene in standard deviation (SD) units. For the GTEx, CMC and ROSMAP data, which are based on RNA sequencing (RNA-Seq), we matched the data sets by Ensembl Gene IDs. For the Braineac data that are based on gene expression microarray, we matched the data sets by gene symbols and removed genes tagged by multiple gene expression probes to ensure a one-to-one match for genes between data sets. The main aim of this study is to quantify the extent to which cis-eQTL data in blood can be used in the SMR analysis21 (see below) to identify genes associated with brain-related phenotypes and disorders. If we had selected the top-associated cis-eQTLs in blood and compared their effects with those in brain, we would likely suffer a form of winner’s curse. To avoid the potential ascertainment bias, we selected the top cis-eQTLs in a reference tissue, i.e., GTEx-muscle (n=361) or CMC (n=467; independent of GTEx), using a stringent Pvalue threshold that is required for the SMR analysis21 (see below), and estimated r b between brain and blood using these SNPs (Supplementary Fig. 3). Although this strategy uses only a quarter of all genes, the estimates of r b should be valid (see below). Note that the estimates of local and distal r g at all SNPs14 would be more informative for other gene-trait association methods such as TWAS22 and MetaXcan23 that use all SNPs in a prediction analysis framework. We chose SMR (URLs) because of one of its features (i.e., the HEIDI test) to filter out associations due to linkage21. We selected the top-associated cis-eQTLs at P eQTL <5×10 −8 for 4257 genes in GTEx-muscle and matched the selected genes with those in the other data sets (the number of matched genes ranged from 1113 to 3841) (Supplementary Table 2, i.e., up to 90%, with the lower numbers matched representing data sets with gene expression data for fewer genes). Note that all the matched genes were expressed in both tissues (i.e., genes which have at least 10 samples with reads per kilobase per million mapped reads (RPKM) > 0.1 and raw read counts >6)13. It should also be noted that our analysis below shows that the test-statistic for the difference in gene expression between tissues was almost independent of the test-statistic for the difference in SNP effect on gene expression between tissues, therefore selecting genes by cis-eQTL P-values would not bias mean gene expression in any specific tissue. We used the Jackknife approach that removes one gene at a time to estimate the sampling variance of ^ rb(Methods) assuming the estimated top cis-eQTL effects for different genes are independent. This assumption was approximately met given the small LD correlations among the 4257 cis-eQTLs and the subtle difference between the mean Jackknife sampling variance and the observed sampling variance in simulation (Supplementary Fig. 4). Results showed that the effects of the top-associated cis-eQTLs were highly correlated between all the brain regions in GTEx after correcting for estimation errors, with a mean ^ rbof 0.94 (s.e. = 0.004; Fig. 1). These estimates are higher than the Spearman correlation estimates reported in a previous study13 because the Spearman correlation does not account for errors in the estimated SNP effects and therefore underestimates the correlation of true effects especially when the sample size is small. The two cerebellum measures (“brain cerebellar hemisphere”and “brain ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04558-1 2NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications cerebellum”) appeared to be outliers. The correlation between “brain cerebellar hemisphere”and “brain cerebellum”was almost perfect (^ rb¼0:99 and s.e. =0.002), but the correlations between the two cerebellum regions and the other regions (mean ^ rb¼0:89 and s.e. =0.006) were significantly smaller than the pairwise correlations between the other regions (mean ^ rb¼0:98 and s.e. =0.003). We performed the same analysis in the Braineac data and observed similar results as above (Supplementary Fig. 5). The estimates of r b between brain and blood in GTEx varied from 0.74 to 0.79 across different brain regions with a mean estimate of 0.77 (s.e. =0.010), similar to the mean estimate of local r g between GTEx-blood and 10 other non-brain GTEx tissues reported in a previous study14. The estimate of r b between CMC (brain) and GTEx-blood was 0.74 (s.e. =0.014), suggesting that the betweensample genetic heterogeneity is small, in line with the strong correlations between CMC and GTEx brain regions (mean ^ rb¼ 0:87 and s.e. =0.010). The estimates of r b from ROSMAP were remarkably similar to those from CMC, providing an important replication of the result. The correlations related to Braineac were notably lower than those related to CMC (Fig. 1), which is likely due to the difference in transcriptomics technology between the two studies (microarray vs. RNA-Seq). It is of note that the results were robust to scale transformation of the eQTL effects (Supplementary Fig. 6), the exclusion of cis-eQTLs in or near the promoter regions (Supplementary Fig. 7), the exclusion of housekeeping genes24,25 (Supplementary Fig. 8), the inclusion of secondary cis-eQTLs identified from a conditional analysis26 (Supplementary Fig. 9), or the adjustment of gene expression data for confounding (e.g., batch effects) predicted from the data (Supplementary Fig. 10). In addition, we selected the topassociated cis-eQTLs at P eQTL <5×10 −8from the CMC data, and found that the estimates of r b among the brain regions and 1 0.96 0.92 0.87 0.83 0.79 0.75 0.7 0.66 0.62 0.58 Brain_Anterior_cingulate_cortex_BA24 Brain_Caudate_basal_ganglia Brain_Cerebellar_Hemisphere Brain_Cerebellum Brain_Cortex Brain_Frontal_Cortex_BA9 Brain_Hippocampus Brain_Hypothalamus Brain_Nucleus_accumbens_basal_ganglia Brain_Putamen_basal_ganglia Whole_Blood CMC ROSMAP Braineac Brain_Anterior_cingulate_cortex_BA24 Brain_Caudate_basal_ganglia Brain_Cerebellar_Hemisphere Brain_Cerebellum Brain_Cortex Brain_Frontal_Cortex_BA9 Brain_Hippocampus Brain_Hypothalamus Brain_Nucleus_accumbens_basal_ganglia Brain_Putamen_basal_ganglia Whole_Blood CMC ROSMAP Braineac 0.974 (0.003) 0.885 (0.007) 0.880 (0.007) 0.990 (0.003) 0.993 (0.002) 0.982 (0.003) 0.980 (0.004) 0.974 (0.003) 0.972 (0.003) 0.763 (0.010) 0.899 (0.010) 0.909 (0.007) 0.767 (0.018) 0.878 (0.007) 0.873 (0.007) 0.968 (0.003) 0.972 (0.003) 0.979 (0.003) 0.970 (0.003) 0.990 (0.002) 0.993 (0.002) 0.766 (0.010) 0.881 (0.009) 0.874 (0.008) 0.762 (0.017) 0.993 (0.002) 0.898 (0.006) 0.891 (0.006) 0.883 (0.007) 0.889 (0.006) 0.874 (0.007) 0.877 (0.007) 0.747 (0.010) 0.818 (0.013) 0.793 (0.012) 0.668 (0.021) 0.903 (0.006) 0.890 (0.006) 0.879 (0.007) 0.887 (0.007) 0.872 (0.007) 0.874 (0.007) 0.743 (0.010) 0.815 (0.013) 0.784 (0.011) 0.661 (0.021) 0.993 (0.002) 0.973 (0.003) 0.971 (0.003) 0.967 (0.003) 0.971 (0.003) 0.764 (0.009) 0.894 (0.009) 0.897 (0.007) 0.757 (0.018) 0.975 (0.003) 0.973 (0.003) 0.973 (0.003) 0.971 (0.003) 0.766 (0.010) 0.904 (0.009) 0.904 (0.007) 0.765 (0.017) 0.978 (0.003) 0.974 (0.003) 0.981 (0.003) 0.792 (0.009) 0.889 (0.010) 0.890 (0.008) 0.773 (0.018) 0.974 (0.003) 0.975 (0.003) 0.782 (0.010) 0.891 (0.010) 0.888 (0.008) 0.778 (0.018) 0.988 (0.002) 0.766 (0.010) 0.886 (0.010) 0.870 (0.008) 0.752 (0.018) 0.775 (0.009) 0.878 (0.010) 0.878 (0.007) 0.769 (0.018) 0.744 (0.014) 0.637 (0.016) 0.584 (0.021) 0.827 (0.012) 0.834 (0.015) 0.688 (0.026) Fig. 1 Estimated correlation of genetic effects of cis-eQTLs between tissues. We estimated r b between brain regions, between brain and blood tissues, and between data sets. The top-associated cis-eQTLs (one for each gene) were selected from GTEx-muscle at P eQTL <5×10 −8. Shown in each cell is the estimate of r b with its standard error given in the parentheses (Methods). In the Braineac data, the eQTLs effect sizes were estimated from gene expression levels averaged across 10 brain regions NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04558-1 ARTICLE NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications 3 between brain and blood in GTEx remained largely unchanged (Supplementary Fig. 11), suggesting that our results are also robust to the ascertainment of the cis-eQTLs. cis-eQTLs with tissue-specific effects. The strong correlation of cis-eQTL effects between brain and blood (Fig. 1) does not preclude eQTLs with detectable difference in effect size between tissues. Of the 1388 cis-eQTLs with P eQTL <5×10 −8in GTExmuscle and available in CMC and GTEx-blood (Supplementary Table 2), 308 (22%) showed significant difference in effect between CMC and GTEx-blood after Bonferroni correction for multiple testing (P difference < 0.05/1388) (Methods). Note that the substantial proportion of eQTLs with significant between-tissue differences in effect does not contradict the large estimate of r b above (Fig. 1) because the power to detect a difference in effect depends on sample size13 (Supplementary Fig. 12). Previous studies have indicated that chromatin state at promoters is largely invariant across diverse cell types whereas enhancers are marked with highly cell-type-specific histone-modification patterns27, that functional variants (predicted by chromatin activity data) in enhancers are less likely to be shared across many tissues compared with those in promoters28, and that cell-type-specific eQTLs are more dispersedly distributed around the transcription start site than eQTLs affected expression in multiple cell types29,30. These results seem to indicate that tissue-specific eQTLs are enriched in distal regulatory elements (i.e., enhancers). To address this hypothesis, we computed the statistics to test for the between-tissue difference in eQTL effect (denoted by T D ) and tested the inflation (or deflation) of mean T D of cis-eQTLs in the functional categories annotated by the Roadmap Epigenomics Mapping Consortium (REMC)31 (Methods). The result showed that although cis-eQTLs are enriched in genomic regions of active chromatin state (e.g., promoters and enhancers) and deflated in inactive regions, the mean T D of cis-eQTLs between CMC and GTEx-blood was almost evenly distributed across all the functional categories with no evidence of inflation in the enhancer regions (Fig. 2). The result remained largely unchanged if we repeated the enrichment analysis based on T D between GTExcerebellum and GTEx-blood (Supplementary Fig. 13). Note that these results do not contradict the observation from a recent study that eQTLs detected in specific tissues in GTEx tend to be most enriched among the variants predicted to be functional in relevant REMC tissues28. There were some examples where the cis-eQTLs with tissue-specific effects in brain and blood were located in enhancers (Supplementary Fig. 14). These examples, however, were rare because only 14 of the 308 eQTLs with P difference < 0.05/1388 were located in enhancers and only 4 of the 14 enhancers appeared to be tissue specific. 0 100 200 300 400 500 Number of eQTLs a 0 1 2 3 4 5 Fold enrichment b TssA Quies ReprPC PromBiv PromP Het ZNFRpts DNase EnhW EnhA TxEn TxWk Tx Prom TssA Quies ReprPC PromBiv PromP Het ZNFRpts DNase EnhW EnhA TxEn TxWk Tx Prom Fig. 2 Enrichment of tissue-specificcis-eQTLs in functional annotations. aThe distribution of cis-eQTLs across 14 functional categories derived from RMEC (Methods). bEstimated enrichment of T D (testing for the difference in cis-eQTL effect between CMC-brain and GTEx-blood) in each functional category (Methods). Error bars represent 95% confidence intervals around the estimates. The black dash line represents fold enrichment of 1. Different colors in a and bcorrespond to 14 functional categories: TssA: active transcription start site, Prom: upstream/downstream TSS promoter, Tx: actively transcribed state, TxWk: weak transcription, TxEn: transcribed and regulatory Prom/Enh, EnhA: active enhancer, EnhW: weak enhancer, DNase: primary DNase, ZNF/ Rpts: state associated with zinc finger protein genes, Het: constitutive heterochromatin, PromP: poised promoter, PromBiv: bivalent regulatory states, ReprPC: repressed Polycomb states, and Quies: a quiescent state 050100150 0 10,000 20,000 30,000 40,000 TD in mean expression level TD in effect size Fig. 3 Correlation of difference in cis-eQTL effect and difference in expression level. Each dot represents one of the 3569 genes between GTEx-cerebellum and GTEx-blood. The 3569 genes were ascertained with at least one cis-eQTL with P eQTL <5×10 −8in GTEx-muscle and expressed in GTEx-cerebellum and GTEx-blood (i.e. genes which have at least 10 samples with RPKM >0.1 and raw read counts >6). In this analysis, we used cis-eQTL effects in SD units and gene expression levels in log 2 (RPKM) units to avoid confounding of the correlation by the mean–variance relationship in gene expression ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04558-1 4NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications In addition, there are a large number of genes showing differences in expression level between tissues11. It is not clear whether these differences are partly driven by the differences in eQTL effect. We sought to address this question by examining the correlation between test-statistic for difference in cis-eQTL effect (in SD units) and test-statistic for difference in mean expression level of the corresponding gene (in log 2 (RPKM) units) between GTEx-cerebellum and GTEx-blood for the 3569 genes each with a cis-eQTL at P eQTL <5×10 −8in GTEx-muscle (Supplementary Table 2). Note that the cis-eQTL effects were re-scaled based on the expression level per gene in SD units so that the correlation was not confounded by the mean–variance relationship in gene expression. That is, if the difference in eQTL effect and that in expression level were both computed in RPKM units, genes with larger differences in mean between tissues are more likely to have differences in variance because of the mean–variance relationship, giving rise to differences in eQTL effect even if the eQTL effects are not different in SD units. We found that the correlation was marginal (r=0.003) (Fig. 3). This is analogous to the observation that there is a large difference in mean height between men and women but the effects of all autosomal SNPs on height in men are almost identical to those in women32,33. However, these results also suggest that an eQTL with identical effect on gene expression in SD units in different tissues could show different effects in RPKM units if the variance of gene expression varies across tissues, which might explain the results from recent studies that genetic variants in or near genes differentially expressed in a particular tissue are enriched for associations with a complex trait34,35. Correlation of cis-mQTL effects between brain and blood. Having shown that cis-eQTL effects are highly correlated between brain and blood, we then turned to estimate the correlation of genetic effects on DNAm between the two tissues by applying the r b method to mQTL data. We analyzed summary-level mQTL data from five studies based on the Illumina HumanMethylation450K array: fetal brain from Hannon et al. (n=166)36, brain cortical region from ROSMAP (n=468)19, frontal cortex region from Jaffe et al. (n=526)37, and peripheral blood from McRae et al. (LBC: n=1366 and BSGS: n=614)38 (Supplementary Table 3). All the mQTL effects are in SD units. We matched the SNPs in common across data sets, selected the top-associated cismQTL at P mQTL <1×10 −10 for 26,840 DNAm probes in the data from Hannon et al. (because only SNPs with P mQTL <1×10 −10 are available in this data set) and matched the selected probes with those in the other data sets (the number of matched probes ranged from 4892 to 6561) (Supplementary Table 4). The correlation of cis-mQTL effects between two brain samples (ROSMAP and Jaffe et al.) was very strong (^ rb¼0:92 and s.e. =0.002), similar to that between two blood samples (^ rb¼0:92 between BSGS and LBC with s.e. =0.003) (Fig. 4a). It is of note that both estimates of r b were smaller than unity, reflecting some degree of heterogeneity between studies. The mean brain–blood r b estimate from two samples was 0.78 (s.e. =0.006) (Fig. 4a), higher than that for cis-eQTLs (mean ^ rb¼0:70 and s.e. =0.015) shown above (Fig. 1). The result remained largely unchanged if the cismQTLs were selected at P mQTL <5×10 −8in the LBC data (Supplementary Fig. 15), again showing the robustness of our results to the choice of reference tissue. In addition, of the 5416 1 0.98 0.95 0.93 0.9 0.88 0.86 0.83 0.81 0.78 0.76 ROSMAP Jaffe et al. BSGS LBC ROSMAP Jaffe et al. BSGS LBC 0.921 (0.002) 0.784 (0.006) 0.773 (0.006) 0.768 (0.006) 0.808 (0.005) 0.923 (0.003) a 0 500 1500 2500 Number of mQTLs b 0.0 1.0 2.0 3.0 Fold enrichment c TssA Quies ReprPC PromBiv PromP Het ZNFRpts DNase EnhW EnhA TxEn TxWk Tx Prom TssA Quies ReprPC PromBiv PromP Het ZNFRpts DNase EnhW EnhA TxEn TxWk Tx Prom Fig. 4 Similarity and difference in cis-mQTL effects between brain and blood. aEstimated r b for cis-mQTLs between brain and blood from four independent data sets. The cis-mQTLs (one for each DNAm probe) were selected at P mQTL <1×10 −10 using data from the Hannon et al. study. Shown in each cell is the estimate of r b with its standard error given in the parentheses (Methods). bThe distribution of cis-mQTLs across 14 functional categories derived from RMEC (Methods). cEstimated enrichment of T D (testing for the difference in cis-mQTL effect between Jaffe-brain and LBC-blood) in each functional category (Methods). Error bars represent 95% confidence intervals around the estimates. The black dash line represents the fold enrichment of 1 NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04558-1 ARTICLE NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications 5 cis-mQTLs, 1847 (34%) showed significantly different effects between brain (Jaffe et al.) and blood (LBC) after correcting for multiple testing (P difference < 0.05/5416). We then tested whether cis-mQTLs in any of the REMC functional categories tend to have higher T D between brain and blood (see above). There were small but significant enrichments of T D in enhancer regions (e.g., transcribed enhancer, active enhancer and weak enhancer) (Fig. 4c), and one of them survived multiple-testing correction (Supplementary Table 5). Meta-analysis of brain eQTL data from correlated samples.We know from the r b analysis above that cis-eQTLs are almost perfectly correlated in different brain regions. We then sought to combine data from the brain regions to increase the power of detecting eQTLs for follow-up analysis (e.g., identification of putative functional genes for brain-related traits and diseases). However, if there is sample overlap between two tissues and the phenotypic correlation is nonzero, the estimation errors of the SNP effects from the two tissues will be correlated. We implemented in the SMR software package (URLs) a summary-databased method, which only requires summary-level data in the cisregions to account for sample overlaps, to meta-analyze cis-eQTL data in correlated samples (MeCS) (Methods). MeCS is very similar to existing meta-analysis approaches such as MTAG39 or the Han et al. method40 that account for sample overlaps. However, there is a small but important distinction. That is, MeCS uses “null”SNPs (e.g., P eQTL > 0.01) to quantify sampling correlation of the estimated SNP effects between two data sets (θ), similar to the strategy used in the latest version of METAL (method unpublished, URLs), whereas MTAG39 uses ^ θestimated by the intercept of bivariate LD score regression41 that relies on the assumption of an infinitesimal model which is invalid in ciseQTL regions42. Han et al.40 suggest the use of the number of overlapping individuals43 or z-statistics to compute ^ θfor summary-data-based analysis. However, a meta-analysis of ciseQTL effects from two tissues requires the correlation of expression level between the tissues (because θ=r p ρwith r p being the correlation of expression level and ρbeing the proportion of sample overlap44) which is not available in summary data, and ^ θestimated by the correlation of z-statistics in the cisregion could be biased by the strong local genetic correlation14. We showed by simulations that ^ θcould be estimated with high accuracy from summary data of the “null SNPs”in cis-region using a simple correlation approach (Supplementary Note 1, Supplementary Figs. 16 and 17), that the MeCS test-statistics were well calibrated under the null hypothesis (Supplementary Fig. 16), and that the MeCS estimates of meta-analysis effect sizes were well estimated under the alternative hypothesis (Supplementary Fig. 17). We compared MeCS to a univariate analysis of the mean expression phenotype across tissues and found that the estimates of effect size and SE from the two approaches were highly consistent (Supplementary Fig. 18). Note that in comparison with the separate analysis in individual tissues, the gain of power for MeCS increased with the decrease of correlation in expression phenotype between tissues, more so for meta-analysis using individuallevel data (Supplementary Fig. 19). We applied MeCS to data from 10 brain regions in GTEx (we referred to the meta-analyzed data as GTEx-brain hereafter). No. of genes with at least one significant cis−eQTL 0 1000 2000 3000 4000 5000 6000 Brain−eMeta ROSMAP CMC GTEx−brain a 56 2041 eQTLGen Brain−eMeta b 33 3030 CAGE Brain−eMeta No. of DNAms with at least one significant cis−mQTL 0 10,000 20,000 30,000 40,000 50,000 BSGS LBC LBC+BSGS Jaffe et al. c 162 34133 d LBC+BSGS Jaffe et al. Fig. 5 Identification of genes and DNAm sites associated with four brain-related traits. Genes (DNAm sites) associated with the brain-related traits were identified by a SMR analysis of GWAS data with eQTL (mQTL) data from brain and blood samples. The four brain-related traits are smoking, IQ, SCZ, and EduYears. a,cshow the number of genes (DNAm sites) with at least one significant SNP at P<5×10 −8in different data sets. b,dshow the number of genes (DNAm sites) associated with traits identified in different data sets. Sample sizes of the brain studies: GTEx-brain (n=~233), CMC (n=467), ROSMAP (n=494), Brain-eMeta (n eff =~1194), and Jaffe et al. (n=526). Sample sizes of the blood studies: CAGE (n=2765), eQTLGen (n=14,115), LBC +BSGS (n=1980) ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04558-1 6NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications There were strong sample overlaps among the ten brain regions (mean overlap =70.4%) and the mean correlation in expression level between pairwise brain regions across all the expressed genes was moderate (mean r p =0.33). The gain of power by the metaanalysis was demonstrated by the observation that the mean χ2 statistic for cis-eQTLs (selected from GTEx-blood at P eQTL <5× 10−8) in GTEx-brain was larger than that in any individual brain region (Supplementary Fig. 20c). The association test-statistic for a SNP can be written as χ2¼1þneff q2 1q2, where n eff is the effective sample size and q2is the variance explained by a SNP45. We therefore can approximately estimate n eff of GTEx-brain assuming constant mean q2across brain regions (Supplementary Note 2). Note that this assumption is justified by the highly consistent estimates of variance of cis-eQTL effects across genes in different brain regions (Supplementary Fig. 21) along with a mean r b estimate of 0.94 between pairwise brain regions for ciseQTL effects in SD units. The estimate of n eff of GTEx-brain was 233, approximately 2.6 times larger than the actual sample size of brain tissue in GTEx (mean n=~89 across 10 brain regions) (Supplementary Fig. 20d). To further increase the power of detecting brain eQTLs, we meta-analyzed GTEx-brain, CMC, and ROSMAP (referred to as Brain-eMeta hereafter). The gain of power is demonstrated by the increased number of genes with at least one cis-eQTL with P eQTL <5×10 −8in Brain-eMeta as compared with that in GTEx-brain, CMC, or ROSMAP (Fig. 5a). Identifying DNAm and genes for brain-related phenotypes. With the Brain-eMeta eQTL data (n eff =1194) obtained from the meta-analysis above, we applied the SMR approach21,46 to test for associations of gene expression levels with four brain-related phenotypes, i.e., ever-smoked (smoking), fluid intelligence score (IQ), years of education (EduYears), and schizophrenia (SCZ). GWAS data were from published meta-analyses for EduYears and SCZ47,48, and from analyses of the full release of the UK Biobank data for smoking and IQ (Methods and Supplementary Table 6). LD data required for the HEIDI test21 were estimated from genotyped/imputed data of the Health and Retirement Study (HRS)49.LDr2from HRS were strongly correlated with those from CMC (Supplementary Fig. 22), consistent with the observation from previous studies26. For power comparison, we included in the SMR analysis an additional set of blood eQTL data from a sample of 14,115 individuals from the eQTLGen Consortium. Only the genes with at least one cis-eQTL at P eQTL < 5×10 −8(one of the basic assumptions of SMR) in both BraineMeta and eQTLGen were included. We further excluded genes in the major histocompatibility complex (MHC) region because of the complexity of this region, leaving 3943 genes for analysis. We identified 61 genes associated with the traits using the brain eQTL data, 41 of which (67.2%) were in common with a larger set of genes (97) identified using the eQTLGen blood eQTL data (Fig. 5b). Despite the heterogeneity between the two eQTL data sets (Brain-eMeta was based on RNA-Seq and eQTLGen was based on microarray), the strong overlap between the two sets of results is consistent with the strong correlation of eQTL effects between brain and blood estimated above. For SCZ, 19 out of the 24 genes identified using brain eQTL data were replicated using blood eQTL data with an additional 27 genes identified only in the blood data because of its larger sample size (Supplementary Fig. 23). We repeated the SMR analysis using blood eQTL data from the Consortium for the Architecture of Gene Expression (CAGE; n=2765)9and observed similar results (Fig. 5b) although the power of CAGE was lower than that of eQTLGen (63 genes identified using CAGE versus 97 genes identified using eQTLGen). We also performed the SMR analysis to detect associations between DNAm sites and the brain-related phenotypes16 using brain mQTL data from Jaffe et al. (n=526) and blood cis-mQTL data from a meta-analysis of LBC and BSGS (n=1980) (Methods). We only included in the analysis DNAm probes with at least one cis-mQTL with P mQTL <5×10 −8in both the brain and blood data sets. We identified 167 DNAm sites associated with the traits (P SMR < 1.8 × 10−6) using the brain mQTL data, 133 of which (79.6%) were in common with the set of 295 DNAm sites identified using the blood mQTL data (Fig. 5d and Supplementary Fig. 24). The brain to blood “replication”rate slightly decreased when we rejected the associations with P HEIDI < 0.05 (Supplementary Fig. 25), likely because of the HEIDI test being over-conservative especially as sample size increases21. These results further demonstrate the feasibility and gain of power of using the cis-genetic effects on gene expression or DNAm estimated in blood to identify putative target genes and regulatory DNA elements for brain-related phenotypes. Discussion We estimated the correlation (^ rb) of genetic effects at the topassociated cis-eQTLs/mQTLs between brain and blood. Because the r b method accounts for estimation errors, ^ rbcan be interpreted as an estimate of correlation of true cis-eQTL effects between brain and blood, as demonstrated by simulations (Supplementary Fig. 2). We applied the method to summary-level eQTL data from GTEx and found that genetic effects on gene expression at the top-associated cis-eQTLs were almost perfectly correlated between different brain regions (mean ^ rb¼0:94 for cis-eQTLs), especially between the non-cerebellar regions (mean ^ rb¼0:98 and s.e. =0.003), in contrast to the modest phenotypic correlation in gene expression levels (mean r p =0.33). It is therefore sensible to run a meta-analysis of the cis-eQTL effects across brain regions to gain power of detecting eQTLs for the whole brain (Supplementary Fig. 19). This can be done even if the brain regions are from different samples. We also found that the cis-eQTL effects were highly correlated between brain and blood in GTEx (mean ^ rb¼0:77 for cis-eQTLs), and the estimate only slightly decreased using data from different samples (mean ^ rb¼0:70). These estimates were significantly different from 1, suggesting there are real genetic differences between tissues. The genetic differences are partly due to cell-type-specific genetic effects regardless whether cell composition covariates have been included in the eQTL analysis or not. This is because adjusting for cell composition only removes the mean differences in gene expression level among cell types rather than cell-type-specific genetic effects. On the other hand, however, the strong betweentissue correlation in cis-eQTL effects does not contradict the result that many genes showed differential expression between brain and blood because the difference in cis-eQTL effect is almost independent of the mean difference in gene expression level (Fig. 3). Our results reinforce that very large sample sizes are needed to generate eQTL data sets in a specific tissue (e.g., blood) to increase the power of detecting cis-eQTLs regardless of the relative expression level of the tissue. Our results also provide some guidelines about the use of discovery–replication paradigm to compare eQTL effects between tissues (i.e., detecting eQTLs in one tissue at a stringent P-value threshold and replicating the effects in another tissue after correcting for multiple tests)13,29. Here, we often saw a low to moderate replication rate even if there is no genetic difference between the tissues. This is because the replication rate is a function of the sample size of the validation set (Supplementary Fig. 12) and the sample sizes of eQTL studies in non-blood tissues are often limited. If we apply the discovery–replication paradigm NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04558-1 ARTICLE NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications 7 to the GTEx data, only ~10.7% of eQTLs discovered in GTExmuscle could be replicated in GTEx-hippocampus (although the estimates from the recent methods50,51 based on the discovery–replication paradigm were much higher) (Supplementary Table 7), which could potentially lead to a wrong conclusion that a large proportion of cis-eQTLs are tissue specific (note that the r b estimate between the two tissues was 0.81). We therefore do not recommend the use of the discovery–replication paradigm to quantify the tissue-specific effects especially in small samples. We applied the SMR and HEIDI methods to identify genes and DNAm sites associated with brain-related phenotypes through pleiotropy using summary data from GWAS and cis-eQTL/ mQTL studies with large sample sizes (n max =453,693 for GWAS, n max =14,115 for eQTL and n max =1980 for mQTL). We identified a number of genes and DNAm sites that showed pleiotropic associations with the phenotypes, consistent with a plausible model that the SNP effects on the phenotypes are mediated by genetic regulation of expression levels of the target genes and/or DNAm levels at the CpG sites. We repeated the analyses using eQTL and mQTL data from brain samples with much smaller sample sizes (n max =1194 for eQTL and n max = 526 for mQTL). Due to the lower power of the data sets, the number of genes or DNAm sites detected in the brain sample was much smaller than that using the blood sample (Fig. 5, Supplementary Figs. 23–25), with at least 50% of genes (DNAm sites) in common between the two sets. These results provide strong justification for the use of blood samples to discover genes related to brain phenotypes and diseases. In practice, we recommend using a blood data set with large sample size for discovery, and an additional data set from brain for replication. This paradigm is certainly applicable to other phenotypes and their related tissues. We conclude with several caveats. First, our estimation of r b is based on genes expressed in both brain and blood (i.e., genes only expressed in one tissue were not included in the estimation). Therefore, the estimate of r b needs to be interpreted with a restriction to genes expressed in both tissues. Although only a quarter (4257) of all genes were selected in our analysis (with at least one cis-eQTL at P eQTL <5×10 −8in GTEx-muscle), up to 90% of those selected genes were expressed in both brain and blood, reflecting the high proportion of all genes expressed in both tissues. Second, we focused our analyses only on cis-eQTLs and cis-mQTLs because trans-eQTLs and trans-mQTLs data were not available in most data sets used in our study. Although most SNP-based heritability for gene expression levels are attributed to cis-eQTLs9,trans-eQTLs may also have an important role in regulating gene expression especially for tissue-specific effects14. The methods developed in this study can be applied to transeQTL/mQTL data with minimal modification. Because the variance explained by individual trans-eQTL/mQTL is small on average9,38, very large sample sizes (e.g., 10,000s) are required to detect trans-eQTLs to be useful for the SMR analysis21. Third, the r b analysis was focused on the correlation at the top-associated cis-eQTLs/mQTLs with relatively large effects (i.e., P<5×10 −8 in a reference tissue) because the SMR test only uses cis-eQTLs/ mQTLs at P<5×10 −8. The estimate of r b was slightly lower for cis-eQTLs/mQTLs selected at a less stringent threshold (Supplementary Fig. 26), consistent with the observation in simulation (Supplementary Fig. 27). However, this does not change our conclusion about the use of the top-associated cis-eQTLs/mQTLs identified in a large blood sample to identify putative target genes for brain-related traits. Last but not least, the MeCS method requires the correlation of errors in the estimated SNP effects between two samples (θ), which is estimated by a simple correlation approach at the null SNPs in the cis-region. This approach, however, is not applicable to eQTL or mQTL summary data that have been ascertained by P-value. It will also be challenging to estimate θif only a small number of cis-SNPs are available in the summary data. We therefore recommend eQTL and mQTL studies to make more cis-SNPs available without ascertainment (e.g., all the cis-SNPs in ±2Mb of a gene or DNAm probe). Despite these caveats, our findings shed light on the genetic architecture underlying the regulation of gene expression across tissues and provide important guidance for studies in the future to identify functional genes for human complex traits. Methods Summary data of cis-eQTL, cis-mQTL, and GWAS. This study is approved by the University of Queensland Human Research Ethics Committee (approval number: 2011001173). All the analyses of eQTL/mQTL data were performed based on summary-level data. A summary description of all the data sets can be found in Supplementary Table 1, Supplementary Table 3, and Supplementary Table 6. All the samples were of European descent and the summary data available to us were derived from individual-level data that passed stringent quantify control (QC) 9,11,18–20,36–38. The SNPs in all eQTL/mQTL data sets were from imputation of the genotyped data to the 1000 Genomes Project (1KGP) reference panels52, and only the SNPs with MAF > 0.01 were included in analyses. The eQTL summary-level data were from six studies, i.e., the Genotype-Tissue Expression (GTEx)11 v6, the CommonMind Consortium (CMC)18, Religious Orders Study and Memory and Aging Project (ROSMAP)19, the Brain eQTL Almanac project (Braineac)20, the Architecture of Gene Expression (CAGE)9, and eQTLGen. In GTEx, ROSMAP, and CMC, gene expression levels were measured by RNA-Seq. Genes in GTEx and ROSMAP were annotated by GENCODE53 v19 and v14, respectively, and genes in CMC were annotated by Ensembl. We accessed the GTEx eQTL summary statistics of ~9.3 million SNPs for ~32,000 genes in 44 tissues (including 10 brain regions) through GTEx portal (URLs). The sample sizes of GTEx tissues ranged from 70 to 361 with an average of 160. We accessed the CMC summary data from Synapse (accession: syn2759792). The CMC eQTL summary statistics (ascertained at FDR < 0.2) of ~1.1 million SNPs for 14,366 genes were derived from individual-level data in dorsolateral prefrontal cortex of 467 subjects, 209 of which were schizophrenia patients. We accessed the ROSMAP eQTL summary statistics of ~6.4 million SNPs for 12,979 genes, which were derived from individual-level data in dorsolateral prefrontal cortex of 494 subjects. We accessed the Braineac eQTL summary statistics of ~6.2 million SNPs for 25,490 genes, which were derived from data in 10 brain regions of 134 subjects free of neurodegenerative disorders20. The gene expression levels in Braineac were measured by Affymetrix Human Exon 1.0 ST Arrays. For blood eQTL data, we used eQTL summary data from CAGE9(38,624 gene expression probes and ~8 million SNPs on 2765 subjects) and eQTLGen (44,556 gene expression probes and ~10 million SNPs on 14,115 subjects). Gene expression levels in CAGE and eQTLGen were measured by Illumina gene expression arrays. We mapped the probes to genes based on the annotations provided by Illumina. The eQTL summary data available in GTEx, CAGE, and eQTLGen were from previous analyses of standardized gene expression levels with mean 0 and variance 1, whereas expression levels in the other data sets (i.e., CMC, ROSMAP, and Braineac) were not standardized. To harmonize the units across data sets, we rescaled the effect size and standard error (SE) of each eQTL in the CMC, ROSMAP, and Braineac based on the z-statistic, allele frequency and sample size using the method described in Zhu et al.21 so that the eQTL effects in all data sets can be interpreted in standard deviation (SD) units. mQTL summary statistics were from five data sets: brain cortical region from ROSMAP study (n ind =468, n probe =420,103, n snp =5 million)19; fetal brain from Hannon et al. (n ind =166, n probe =26,840, n snp =0.3 million)36; frontal cortex region from Jaffe et al. (n ind =526, n probe =138,917, n snp =1.5 million)37; and peripheral blood from McRae et al.38 (Lothian Birth Cohorts54 (LBC): n ind =1366 and Brisbane Systems Genetics Study55 (BSGS): n ind =614). DNAm levels in all these five studies were based on the Illumina HumanMethylation450K array. We performed a meta-analysis of LBC and BSGS, resulting in 397,621 DNAm probes and ~7.7 million SNPs. The DNAm levels of all the five studies were not standardized. We computed the effect size and SE of each mQTL from their zstatistics using the method described in Zhu et al.21 We included in the analysis four brain-related complex traits, i.e., ever-smoked (smoking), fluid intelligence score (IQ), years of education (EduYears), and schizophrenia (SCZ). GWAS summary statistics for EduYears (n=293,723) and SCZ (36,989 cases and 113,075 controls) were from the latest meta-analyses47,48, and summary data for smoking (n=453,693) and IQ (n=146,819) were from GWAS analyses of the latest release of the UK Biobank (UKB) data56. Quality control and imputation of the UKB data have been detailed elsewhere56. We used 456,426 individuals of European descent and 7,288,503 common SNPs (MAF > 0.01) imputed from the Haplotype Reference Consortium (HRC)57 reference panel in the analysis. IQ was measured by 13 fluid intelligence questions and detailed description of the measurement can be found at the UKB website (URLs). We adjusted IQ (n=146,819) by age and sex, and standardized the adjusted phenotype by rank-based inverse-normal transformation. The GWAS analyses were ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04558-1 8NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications performed in BOLT-LMM58 using all 7.3 million SNPs with a subset of 0.7 million SNPs in common with HapMap359 used to control for population structure and polygenic effects. We used self-reported “ever-smoked”as a dichotomous phenotype for smoking (208,988 cases and 244,705 controls). We analyzed the data in BOLT-LMM based a linear model with age and sex fitted as covariates, and transformed the effect size of each SNP on the observed 0–1 scale to odds ratio (OR) using LMOR60 (URLs). Correlation of cis-eQTL effects between tissues. Let ^ bbe the estimated effect at the top-associated cis-eQTL for a gene (i.e., one SNP per gene). We can model ^ bas ^ b¼bþeð1Þ where bis the true effect and eis the estimation error. We assume that band eare random variables when interrogated across genes, i.e., bN0;var bðÞðÞand eN0;var eðÞðÞ. The covariance of the estimated cis-eQTL effects between tissues i and jacross genes can be partitioned into the covariance of true cis-eQTL effects and the covariance of estimation errors (if there is a sample overlap), i.e., covð ^ bi; ^ bjÞ¼covðbi;bjÞþcovðei;ejÞ¼covðbi;bjÞþreffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi varðeiÞvarðejÞ qð2Þ where var ei ðÞ and varðejÞare the variance of the estimation errors across genes in tissues iand j, respectively, and r e is the correlation of estimation errors across genes between two tissues, i.e., re¼corðei;ejÞ. We know from Bulik-Sullivan et al.41 and Zhu et al.44 that rerpρ, where ρ¼Ns ffiffiffiffiffiffiffi NiNj pmeasures the sample overlap with N i and N j being the sample sizes in tissues iand j, respectively, and N s being the number of overlapping individuals, and r p is the correlation of gene expression levels between two tissues in the overlapping sample. If i=j, then r e =1 and var bi ðÞ¼varð ^ biÞvar ei ðÞ, where var(b i ) is the variance of true cis-eQTL effects across genes in tissue i. We therefore can estimate the correlation of true cis-eQTL effect sizes across genes between tissues iand jas ^ rb¼c cov bi;bj  ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c var bi ðÞ c varðbjÞ q¼c cov ^ bi; ^ bj  ^ reffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c var ei ðÞ c varðejÞ q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c var ^ bi  c var ei ðÞ hi c var ^ bj  c var ej hi rð3Þ where c varð ^ biÞand c varð ^ bjÞ(i.e., the estimates of varð ^ biÞand varð ^ bjÞ) are the observed sample variances of ^ biand ^ bj, respectively, in a set of genes, and c covð ^ bi; ^ bjÞis the observed sample covariance between ^ biand ^ bjin the set of genes. However, c var ei ðÞ,c varðejÞand ^ reare not observable. We know that SE2of ^ bof a SNP is an estimate of the variance of eover repeated experiments for a gene. We therefore can approximate c varðeÞby the average of SE2across genes (one SNP per gene). We also know from Eq. (2) that if b i =b j =0, covð ^ bi; ^ bjÞ¼reffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi var ei ðÞvarðejÞ q. Hence, ^ re¼b cov ^ bi; ^ bj ðÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b var ei ðÞ b varðejÞ p¼b cov ^ bi; ^ bj ðÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi b var ^ bi ðÞ b var ^ bj ðÞ p¼ c corð ^ bi; ^ bjÞfor null SNPs, where c corð ^ bi; ^ bjÞis the observed sample correlation between ^ biand ^ bjin the set of genes. In practice, we computed ^ refor each gene using “null”SNPs (P eQTL > 0.01) in the cis-region by a simple correlation approach and took the average across genes. The sampling variance of ^ rbover repeated experiments can be computed via Jackknife approach leaving one gene out at a time. ^ V^ rb ðÞ Jackknife¼m1 mX t ^ rbðtÞ^ rbð:Þ hi 2 ð4Þ where ^ rbðtÞis the estimate with the t-th gene left out and ^ rbð:Þ¼1 mP t ^ rbðtÞ. The method is derived based on eQTL data but can be applied to data from genetic studies of different types of molecular phenotypes (e.g., DNAm and histone modification). Enrichment of tissue-specific eQTLs in functional categories. We used chromatin state data from 23 blood samples (T-cell, B-cell, and hematopoietic stem cells) and 10 brain samples generated by the NIH Roadmap Epigenomics Mapping Consortium (REMC)31. There were 25 chromatin states predicted by ChromHMM61 based on the imputed data of 12 histone-modification marks31.We classified the 25 chromatin states into 14 main functional categories by combining functionally relevant annotations62. We tested the difference in eQTL effect for a gene between two tissues (iand j) using the method below. Let ^ d¼ ^ bi ^ bjð5Þ The sampling variance of ^ dover repeated experiments can be written as V ^ d ¼V ^ bi  þV ^ bj  2θffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V ^ bi  V ^ bj  rð6Þ where ^ biand ^ bjare the estimated effect sizes of the top-associated cis-eQTL for a gene in two tissues iand j,Vð ^ biÞand Vð ^ bjÞare the sampling variances of ^ biand ^ bj, respectively, over repeated experiments, and θis sampling correlation between ^ bi and ^ bjfor the gene over repeated experiments. In practice, ^ Vð ^ biÞand ^ Vð ^ bjÞcan be estimated by SE2of ^ biand ^ bj, and ^ θcan be approximated by the sample correlation between ^ biand ^ bjacross the “null”SNPs (e.g., P eQTL > 0.01) in the cis-region for the gene. The significance of ^ dcan therefore be assessed by a Wald test, i.e., TD¼ ^ d2 var ^ d ðÞχ2 1. To test the enrichment of T D statistics in functional annotations, we allocated the cis-eQTLs to the 14 functional categories described above by physical position, and calculated the mean T D of each category. We assessed the enrichment using the inflation factor λ=TDðiÞ TD , where TDðiÞis the mean T D of the cis-eQTLs in a category i, and TDis the mean T D of all the cis-eQTLs. We then used the Jackknife approach (leaving one gene out at a time) described above to compute the variability of λ. Note that although we described the enrichment test method above based on ciseQTLs, the method can be applied to data from genetic studies of different types of molecular phenotypes (e.g., DNAm and histone modification). Meta-analysis of cis-eQTL data from correlated samples. We know from Eq. (1) that the estimated effect of a cis-eQTL for a gene can be partitioned into two components, i.e., the true effect size (b) and the estimation error (e). For multiple tissues, the joint distribution of the estimates can be written as ^ bN1b;SðÞ ð7Þ where ^ b¼½ ^ b1; ^ b2;¼; ^ bt,Sis the sampling (co)variance matrix of ^ bover repeated experiments with Sij ¼Cð ^ bi; ^ bjÞ.Sij ¼θijSiSjwhen i≠j, where θ ij is sampling correlation between ^ biand ^ bjfor the gene over repeated experiments. S2 iand S2 jare the sampling variance of b biand b bj, respectively, over repeated experiments. If i=j, then θij ¼1 and Sij ¼S2 i. In practice, ^ θij can be approximated by the sample correlation of the estimated SNP effects between a pair of tissues across the “null” SNPs (e.g., P eQTL > 0.01) in the cis-region for each gene. Similar to the summarydata-based meta-analysis methods that account for correlated estimation errors39,40,63, we can estimate combined effect as ^ b¼1T^ S11  11T^ S1^ bð8Þ ^ V ^ b ¼1 1T^ S11ð9Þ The significance of ^ bcan be assessed by a Wald test, i.e., ^ b2 ^ V ^ b ðÞχ2 1. URLs. For MeCS, see http://cnsgenomics.com/software/smr/#MeCS. For SMR, see http://cnsgenomics.com/software/smr. For LMOR, see http://cnsgenomics.com/ shiny/LMOR/. For UK Biobank, see http://biobank.ctsu.ox.ac.uk/. For METAL, see https://genome.sph.umich.edu/wiki/METAL. For GTEx Portal, see http://www. gtexportal.org/. For CMC data, see https://www.synapse.org/CMC. For Braineac data, see http://www.braineac.org/. Data availability. Brain-eMeta eQTL summary data are available at http:// cnsgenomics.com/software/smr/#Download. The eQTLGen summary data are available through application to the eQTLGen consortium. All the other data sets used in this study are from the public domain. The software tools are available at the URLs above. Received: 22 February 2018 Accepted: 10 May 2018 References 1. Welter, D. et al. The NHGRI GWAS Catalog, a curated resource of SNP-trait associations. Nucleic Acids Res. 42, D1001–D1006 (2013). 2. Visscher, P. M., Brown, M. A., McCarthy, M. I. & Yang, J. Five years of GWAS discovery. Am. J. Hum. Genet. 90,7–24 (2012). 3. Visscher, P. M. et al. 10 years of GWAS discovery: biology, function, and translation. Am. J. Hum. Genet. 101,5–22 (2017). 4. Roussos, P. et al. A role for noncoding variation in schizophrenia. Cell Rep. 9, 1417–1429 (2014). 5. Farh, K. K.-H. et al. Genetic and epigenetic fine mapping of causal autoimmune disease variants. Nature 518, 337–343 (2015). NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-04558-1 ARTICLE NATURE COMMUNICATIONS | (2018) 9:2282 |DOI: 10.1038/s41467-018-04558-1 |www.nature.com/naturecommunications 9