A Statistical Framework for Joint eQTL Analysis in Multiple Tissues

Summary

Mapping expression Quantitative Trait Loci (eQTLs) represents a powerful and widely adopted approach to identifying putative regulatory variants and linking them to specific genes. Up to now eQTL studies have been conducted in a relatively narrow range of tissues or cell types. However, understanding the biology of organismal phenotypes will involve understanding regulation in multiple tissues, and ongoing studies are collecting eQTL data in dozens of cell types. Here we present a statistical framework for powerfully detecting eQTLs in multiple tissues or cell types (or, more generally, multiple subgroups). The framework explicitly models the potential for each eQTL to be active in some tissues and inactive in others. By modeling the sharing of active eQTLs among tissues, this framework increases power to detect eQTLs that are present in more than one tissue compared with “tissue-by-tissue” analyses that examine each tissue separately. Conversely, by modeling the inactivity of eQTLs in some tissues, the framework allows the proportion of eQTLs shared across different tissues to be formally estimated as parameters of a model, addressing the difficulties of accounting for incomplete power when comparing overlaps of eQTLs identified by tissue-by-tissue analyses. Applying our framework to re-analyze data from transformed B cells, T cells, and fibroblasts, we find that it substantially increases power compared with tissue-by-tissue analysis, identifying 63% more genes with eQTLs (at FDR = 0.05). Further, the results suggest that, in contrast to previous analyses of the same data, the majority of eQTLs detectable in these data are shared among all three tissues.

Introduction

Regulatory variation plays an essential role in the genetics of disease and other phenotypes as well as in evolutionary change [1]–[3]. However, in sharp contrast to nonsynonymous variants in the human genome, which can now be identified with great accuracy, it remains extremely difficult to know which variants in the genome may impact gene regulation in any given tissue or cell type. [Henceforth we use “tissue” for brevity, but everything applies equally to cell types.] Expression QTL mapping (e.g. [4]–[6] represents a powerful approach for bridging this gap, by allowing regulatory variants to be identified, and linked to specific genes. Indeed, numerous studies (e.g. [7], [8]) have shown highly significant overlaps between eQTLs and SNPs associated with organismal-level phenotypes in genome-wide association studies (GWAS), suggesting that a large fraction of GWAS associations may be due to variants that affect gene expression.

Ultimately, understanding the biology of organismal phenotypes, such as diseases, is likely to require understanding regulatory variation in many different tissues ([9], [10]). For example, if regulatory variants differ across tissues, then, in understanding GWAS hits, and using them to understand the biology of disease, we would like to know which variants are affecting which tissues. At a more fundamental level, identifying differential genetic regulation in different tissues could yield insights into the basic biological processes that influence tissue differentiation. To date, eQTL studies have been performed in a relatively narrow range of tissue types. However, this is changing quickly: for example, the NIH “Genotype-Tissue Expression” (GTEx) project aims to collect expression and genotype data in 30 tissues across 900 individuals. Motivated by this, here we describe and illustrate a statistical framework for mapping eQTLs in expression data on multiple tissues.

While statistical methods for identifying eQTLs in a single tissue or cell type are now relatively mature (e.g. [11]) current analytic tools are limited in their ability to fully exploit the richness of data across multiple tissues. In particular, available methods fall short in their ability to jointly analyze data on all tissues to maximize power, while simultaneously allowing for differences among eQTLs present in each tissue. Indeed relatively few papers have considered the problem. The simplest approach (e.g. [12], [13]) is to analyze data on each tissue separately (“tissue-by-tissue” analysis), and then to examine overlap of results among tissues. However, this fails to leverage commonalities among tissues to improve power to detect shared eQTLs. Furthermore, although examining overlap of eQTLs among tissues may appear a natural approach to examining heterogeneity, in practice interpretation of results is complicated by the difficulty of accounting for incomplete power. Both [13] and [14] provide approaches to address this, but only for pairwise comparisons of tissues.

Compared with tissue-by-tissue analysis, joint analysis of multiple tissues has the potential to increase power to identify eQTLs that have similar effects across tissues. Both [15] and [16] conduct such joint analyses – the first using ANOVA, and the second using a weighted -score meta-analysis – and [16] confirm that their joint analysis has greater power than tissue-by-tissue analysis. The ANOVA and -score methods each have different advantages. The ANOVA framework has the advantage that, by including interaction terms, it can be used to investigate heterogeneity in eQTL effects among tissues. Gerrits et al. ([15]) use this to identify eQTLs that show significant heterogeneity, and then classify these eQTLs, post-hoc, into different types based on estimated effect sizes. The weighted -score method has the advantage that, unlike ANOVA, it allows for different variances of expression levels in different tissues (which are likely to occur in practice). However, it does not so easily allow for investigation of heterogeneity; Fu et al. ([16]) hence assess heterogeneity for pairs of tissues by using a resampling-based procedure to assess the significance of observed differences in scores. Other papers, including [17] and [18], also show that joint analyses provide more power. Our work goes beyond these papers in its modeling of heterogeneity, and in its use of a hierarchical model to borrow information across genes to estimate weights associated with different types of heterogeneity.

Here we introduce a statistical framework for the joint analysis of eQTLs among multiple tissue types, that combines advantages of some of the methods above, as well as introducing some new ones. In brief, our framework integrates recently-developed GWAS meta-analysis methods that allow for heterogeneity of effects among groups [19]–[23], into a hierarchical model (e.g. [24], [25]) that combines information across genes to estimate the relative frequency of patterns of eQTL sharing among tissues. Like ANOVA, our approach allows investigation of heterogeneity among several tissues, not just pairs of tissues. However, in contrast to ANOVA, our framework allows for different variances in different tissues. Moreover, unlike any of the methods described above, our framework explicitly models the fact that some tissues may share eQTLs more than others, and estimates these patterns of sharing from the data (a similar idea was applied to ChIP-Seq data by [26]). Our methods also allow for intra-individual correlations when samples are obtained from a common set of individuals. While we focus here on comparing and combining information across different tissue types, our framework could be applied equally to comparing and combining across other units, e.g. different experimental platforms, multiple datasets on the same tissue types, or data on individuals from different populations.

The remainder of the paper is as follows. After providing a brief overview of our framework, we use simulations to illustrate its power compared to other methods, and then apply it to map eQTLs, and assess heterogeneity among tissues, using data from Fibroblasts, LCLs and T-cells ([12]). Consistent with results from [16], we show that our joint analysis framework provides a large gain in power compared with a tissue-by-tissue analysis. Furthermore, compared with previous analyses of these data, we find a much higher rate of tissue-consistent eQTLs.

Results

Methods Overview

Consider mapping eQTLs in tissues. In our applications here the expression data are from micro-arrays, and so we assume a normal model for the expression levels, suitably-transformed. (These methods can also be applied to RNA-seq data after suitable transformation; see Discussion). That is, in each tissue, , we model the potential association between a candidate SNP and a target gene by a linear regression:

where denotes the observed expression level of the target gene in tissue for the individual, the mean expression level of this gene in tissue , the effect of a candidate SNP on this gene expression in tissue , the genotype of the individual at the SNP (coded as 0,1 or 2 copies of a reference allele) and the residual error for tissue and individual . Note that the subscript on residual variance indicates that we allow the residual variance to be different in each tissue. In addition, when tissues are sampled from the same set of individuals, we allow that the residual errors may be correlated (with the correlation matrix to be estimated from the data).

The primary questions of interest are whether the SNP is an eQTL in any tissue, and, if so, in which tissues. To address these questions we use the idea of a “configuration” from [21], [23]. A configuration is a binary vector where indicates whether the SNP is an eQTL in tissue . If then we say the eQTL is “active” in tissue . The “global null hypothesis”, , that the SNP is not an eQTL in any tissue, is therefore . Every other possible value of can be thought of as representing a particular alternative hypothesis. For example, represents the alternative hypothesis that the SNP is an eQTL in all tissues, and represents the alternative hypothesis that the SNP is an eQTL in just the first tissue.

Our aim is to perform inference for . A natural approach is to specify a probability model, , being the probability of obtaining the observed data if the true configuration were , and then perform likelihood-based inference for . The support in the data for each possible value of , relative to the null , is quantified by the likelihood ratio, or Bayes Factor (BF, [27]):

Specifying these likelihoods requires assumptions about , the distribution of the effect sizes for each possible configuration (as well as less crucial assumptions about nuisance parameters such as and ). Of course, if then by definition, but for the tissues where various assumptions are possible – for example, one could assume that the effect is the same in all these tissues, or allow it to vary among tissues. Here we use a flexible family of distributions, (see Methods), where the hyper-parameters can be varied to control both the typical effect size, and the heterogeneity of effects across tissues (see below).

The value of measures the support in the data for one specific alternative configuration , compared against the null hypothesis . To account for the fact that there are many possible alternatives, the overall strength of evidence against at the candidate SNP is obtained by “Bayesian Model Averaging” (BMA), which involves averaging over the possible alternative configurations , weighting each by its prior probability, :

Further, under an assumption of at most one eQTL per gene, the overall evidence against for the entire gene (i.e. that the gene contains no eQTL in any tissue) is given by averaging across all candidate SNPs [28]. In either case, at either the SNP or gene level, large values of constitute strong evidence against . has a direct Bayesian interpretation as the strength of the evidence against , but here we also use it as a frequentist test statistic ([28], [29]), assessing significance by permutation or simulation. The latter has the advantage that and obtained in this way are “valid” even if not all the prior assumptions are exactly correct.

Note that depends on the choice of , and the power of as a test statistic is expected to depend on how well this choice of these hyper-parameters captures the range of alternative scenarios present in the data. Here we make use of three different choices:

• A “data-driven” choice, where the hyper-parameters are estimated from the data using a hierarchical model (HM, [30]) that combines information across all genes. We use to denote this choice.

• A “default” choice, which chooses to cover a wide range of different possible alternative configurations, and is set to allow modest heterogeneity. We use to denote this choice.

• A “lite” choice, which puts weight only on the most extreme configurations (where the eQTL is active in only one tissue, or in all tissues), but compensates by setting to allow for more heterogeneity. We use to denote this choice.

Each of these choices has something to recommend it. The first, being data driven, is the most attractive in principle, but also the most complex to implement. The default choice is simpler to implement, and is included partly to demonstrate that one does not have to get the hyper-parameter values exactly “right” for to be a powerful test statistic. Finally, has the advantage that it is easily applied to large numbers of tissues; neither of the other methods scales well, either computationally or statistically, with the number of tissues, because the number of terms in the sum in equation (3) is .

When there is strong evidence against , the Bayes Factors can also be used to assess which alternative configurations are consistent with the data. Specifically the posterior probability on each configuration is:

and the posterior probability that the SNP is an eQTL in tissue is obtained by summing the probabilities over configurations in which :

The second of these is particularly helpful when the data are informative for an eQTL in tissue , but ambiguous in other tissues: in such a case the probability will be close to 1, even though the “true” configuration will be uncertain (so none of the probabilities (4) will be close to 1). Because both (4) and (5) are sensitive to choice of hyper-parameters, we compute them using (where the hyper-parameters are estimated from the data).

Further details of methods used are provided in the Methods section.

Simulations

Power to detect eQTLs

We begin by comparing the ability of different methods to reject the global null hypothesis ; i.e. to detect eQTLs that occur in any tissues. We expect that a tissue-by-tissue analysis, which analyzes each tissue separately, will perform well for detecting eQTLs that are present in a single tissue. Conversely, we expect joint analysis of all tissues to perform well for detecting eQTLs that are present in all tissues. Our Bayesian model averaging (BMA) approach attempts, by averaging over different possible eQTL configurations, to combine the advantages of both types of analysis, and thus aims to perform well across all scenarios.

To assess this we performed a series of simulations, with five tissues measured in 100 individuals (and no intra-individual correlations). Each simulation consisted of 2,000 gene-SNP pairs (i.e. one candidate SNP per gene), half of which were “null” (i.e. the SNP was not an eQTL in any tissue), and the other half following an alternative hypothesis where the SNP was an eQTL in exactly tissues, with varying from 1 to 5. Thus, for example, the simulations with assess power to detect eQTLs that are active in just one tissue, whereas the simulations with assess power to detect eQTLs that are active in all five tissues. When simulating eQTLs that are active in multiple tissues we assumed their effects to be similar, but not identical, across tissues (see Methods). We applied four analysis methods to these data: 1) , being our Bayesian Model Averaging approach with default weights described above; 2) , being the computationally-scalable version of BMA described above; 3) ANOVA/linear regression (ANOVA/LR) (c.f. [15] and see Methods), which jointly analyzes all tissues in a regression model, and compares the general alternative model (which allows a different genetic effect in each tissue) with the null model (no effect in all tissues); and 4) a “tissue-by-tissue” analysis (c.f. [12]), where we use linear regression to test for an eQTL separately in each tissue, and take the minimum across tissues as a test statistic. For simplicity we defer consideration of the more sophisticated of our approaches, , to slightly more complex simulations described later. Each of these methods produces a test statistic for each SNP-gene pair, testing the global null hypothesis . For each test statistic, we found the threshold that yielded a False Discovery Rate of 0.05 (based on the known null/alternative status of each SNP-gene pair), and assessed the effectiveness of each method by the number of discoveries it made at that FDR.

The results of these comparisons are shown in Figure 1A. As expected, for eQTLs that occur in just one tissue, the tissue-by-tissue analysis performs best. However, it is only slightly more effective than the joint analysis approaches in this setting. Conversely, the joint analysis approaches outperform the tissue-by-tissue analysis for eQTLs that occur in more than one tissue, with the gains becoming larger as the number of tissues sharing the eQTL increases. The BMA analyses generally perform similarly to one another, and outperform ANOVA/LR. This is presumably because our simulations involved eQTLs that have similar effects in each tissue, and our prior distribution explicitly up-weights eQTLs with this feature.

This first set of simulations assumed error variances to be equal among tissues. This assumption is made by ANOVA/LR, but not by the other methods, and is likely often to be violated in practice. To assess the effects of this we repeated the simulations, but with error variances differing among tissues. The results (Figure 1B) confirm that, relative to other methods, ANOVA/LR performs less well when error variances vary among tissues.

To assess performance in larger numbers of tissues we repeated the simulations above, but with 20 tissues (so ). For this many tissues computing involves averaging over all possible eQTL configurations, which is computationally inconvenient, so we omitted from this comparison. The results (Figure 1C) show that performs similarly to the tissue-by-tissue analysis for eQTLs that occur in just one or two tissues, and outperforms it substantially for eQTLs occurring in many tissues. As expected, ANOVA/LR outperforms tissue-by-tissue analysis for eQTLs occurring in many tissues, but is noticeably less effective for eQTLs occurring in only one tissue, and performs consistently less well than .

In summary, these simulations illustrate the benefits of Bayesian Model Averaging as a general strategy for producing powerful test statistics: by explicitly averaging over a range of alternative models, these test statistics are able to achieve good power to detect a wide range of different types of eQTL.

Identifying eQTLs in particular tissues: borrowing information among tissues

Next we consider the benefits of jointly analyzing multiple tissues, even when the main goal is to identify eQTLs in a particular tissue of interest. For intuition, suppose for a moment that every eQTL is shared among all tissues. Then, from the simulation results above, we know that a joint analysis will identify more eQTLs overall, and hence more eQTLs in the tissue of interest. Of course, not all eQTLs are shared among all tissues, but some are, and some tissues may share eQTLs more than others. To allow for this, our hierarchical model attempts to infer the extent of such sharing (by estimating the configuration weights ), and exploits any sharing that does occur to increase power to detect eQTLs in each tissue. By estimating sharing from the data, rather than assuming that all tissues share equally with one another (as do the simpler test statistics and used above), we expect to make more effective use of sharing in the data to further improve power to identify eQTLs.

To illustrate this, we simulated eQTL data for five tissues. Some eQTLs were shared by all tissues, some were specific to each tissue, some were shared by Tissues 1 and 2 only, and some were shared by Tissues 3, 4 and 5. To show how the benefits of sharing can change with sample size, we simulated 60 samples for Tissue 1, and 100 samples for the others. This mimics a setting where Tissue 1 is harder to obtain than the other tissues, with Tissue 2 being the best proxy for Tissue 1.

We applied our Bayesian methods and a tissue-by-tissue analysis to these data, and assessed their ability to identify eQTLs in each tissue. For the tissue-by-tissue analysis the test statistic in each tissue is simply the linear regression in that tissue. For our Bayesian methods, the test statistic in each tissue is the posterior probability of the SNP being an active eQTL in that tissue (5). Note that this posterior probability is computed from joint analysis of all tissues, and takes account of sharing of eQTLs among tissues. For example, consider a SNP showing modest association with expression in Tissue 1. If this SNP also shows strong association in the other tissues, then it will be assigned a higher probability of being an active eQTL in Tissue 1 than it would if it showed no association in the other tissues. For each method, separately in each tissue, we identified the threshold of the test statistic value that yields a FDR of 0.05 in that tissue, based on the true active/inactive status of each SNP in that tissue (known since this is simulated data). (A SNP that is an eQTL in some tissues but not others counts as a “false discovery” if it is called as an eQTL in a tissue where it is inactive.) For the Bayesian methods we obtained results both using “default” weights on configurations (), and using weights estimated from the data by the hierarchical model (). The latter is expected to be more effective as it should learn, for example, that Tissue 1 shares more eQTLs with Tissue 2 than with other tissues.

The results (Figure 2) show that, for all tissues, our joint analyses outperform the tissue-by-tissue analysis. Further, outperforms , demonstrating the benefits of learning patterns of sharing from the data. Finally, the gain of is greater for Tissue 1 than for Tissue 2, illustrating that benefits of sharing information are greater for tissues with small sample sizes.

Furthermore, using the hierarchical model which pools all genes together, we can estimate the configuration proportions, i.e. . In the setting described above, we simulated one thousand eQTLs in each of 8 different configurations, as well as one thousand genes with no eQTLs. Averaged over replicates, the proportions are estimated to be in for each of the 8 active configurations (negligible differences between replicates). These estimates are fairly accurate knowing that the true proportion is for each configuration.

Analysis of eQTL Data in Three Cell Types from Dimas et al.

We now analyze data from [12], consisting of gene expression levels measured in fibroblasts, LCLs and T-cells from 75 unrelated individuals genotyped at approximately 400,000 SNPs. The data were pre-processed similarly to the original publication, as described in the Methods section. Throughout we focus on testing SNPs that lie within 1 Mb of the transcription start site of each gene (the “cis candidate region”), and on a subset of 5012 genes robustly expressed in all three cell-types.

Gain in power from joint analysis

First we assess the gain in power from mapping eQTLs in all three cell types jointly, using , compared with a “tissue-by-tissue” analysis similar to that in [12]. For each method we compute a test statistic for each gene, combining information across SNPs, to assess the overall support for any eQTL in that gene in any tissue. For our Bayesian approach the test statistic is the average value of over all SNPs in the cis candidate region; for the tissue-by-tissue analysis the test statistic is the minimum from linear regressions performed separately in each tissue for each SNP in the cis candidate region. We translate each of these test statistics into a for each gene by comparing observed values with simulated values obtained under (by permutation of individual labels). Finally we computed, for each method, the number of genes identified as having an eQTL at an FDR of 0.05 (using the qvalue package [31]).

Joint mapping, via , substantially increased power to identify eQTLs compared with tissue-by-tissue analysis. For example, identified 1022 eQTLs at FDR = 0.05, which is more than the 627 eQTLs identified by the tissue-by-tissue analysis at the same FDR (Figure 3A, 3B, and 3C). Further, the vast majority of eQTLs identified by the tissue-by-tissue analysis () are also detected by (Figure 3C).

In many cases, the eQTLs detected by but not by the tissue-by-tissue analysis have modest effects that are consistent across tissues. Because their effects are modest in each tissue, they fail to reach the threshold for statistical significance in any single tissue, and so the tissue-by-tissue analysis misses them. But because their effects are consistent across tissues, the joint analysis is able to detect them. Figure 4 shows an example illustrating this (gene ASCC1, Ensembl id ENSG00000138303, with SNP rs1678614). The PC-corrected phenotypes already indicate that this gene-SNP pair looks like a consistent eQTL (Figure 4A), and its effect size estimates are highly concordant across tissues (Figure 4B). However, as indicated by the on the forest plot, this eQTL is not called by the tissue-by-tissue analysis in any tissue (all the exceed ). In contrast, the joint analysis pools information across tissues to conclude that there is strong evidence for an eQTL (), and that it is likely an eQTL in all three tissues (probability 1 assigned to the consistent configuration , conditional on it being an eQTL).

Many eQTLs are consistent among tissues

The original analyses of these data concluded that to of eQTLs operate in a cell-type specific manner ([12]). These results were obtained by mapping eQTLs separately in each tissue, and then examining which of the eQTLs identified in each tissue also showed some signal (e.g. at a relaxed significance threshold of ), in another tissue. However, as noted by [14], due to incomplete power, eQTLs that are actually shared between tissues may appear “tissue-specific” in this type of analysis. Our hierarchical model has the potential to help overcome this difficulty by estimating the proportion of eQTLs that follow each configuration type as a parameter of the model, combining information across all genes, and without setting specific significance thresholds (thus sidestepping the problems of incomplete power).

Applying the hierarchical model to these data produced an estimate of just of eQTLs being specific to a single tissue, with an estimated of eQTLs being common to all three tissues (95% CI = 84%–93%; Table 1). Among eQTLs shared between just two tissues, many more are shared between LCLs and T-cells (Figure S1 shows such an example), than between these cell types and fibroblasts. This is consistent with results from [12], and perhaps expected since LCLs and T-cells are more similar to one another than to fibroblasts.

We obtained qualitatively similar patterns when we varied some of the assumptions in the hierarchical model - specifically, whether or not we allow for intra-individual correlations, whether or not we assume at most one eQTL per gene, whether or not we remove PCs to account for confounders, and whether or not we analyze all genes or only those genes robustly expressed in all tissues (Text S1). Nonetheless, we caution against putting too much weight on any particular number to quantify tissue specificity, not least because the definition of a tissue-specific eQTL is somewhat delicate: for example, it is unclear how to classify a SNP that is very strong eQTL in one tissue, and much weaker in the others. Further, our estimates necessarily reflect patterns of sharing only for moderately strong eQTLs, strong enough to be detected in the modest sample sizes available here: patterns of sharing could be different among weaker eQTLs. Nonetheless, these results do suggest that there is substantial sharing of eQTLs among these three tissue types, considerably higher than the original analysis suggested.

To illustrate the potential pitfalls of investigating heterogeneity in a tissue-by-tissue analyses, we also ran a tissue-by-tissue analysis on these data. Specifically, we called eQTLs separately in each tissue (at an FDR of 0.05), and then examined the overlap in the genes identified in each tissue. Using this procedure, in strong contrast with results from the joint analysis, of eQTLs are called in only one tissue, with fewer than 15% called in all three tissues (Table 1). Qualitatively similar results are obtained for different FDR thresholds. However, these results cannot be taken as reliable indications of tissue specificity, because the procedure fails to take account of incomplete power to detect eQTLs at any given threshold, and therefore tends to over-estimate tissue specificity. Figure 5 shows an eQTL that illustrates this behavior (gene CHPT1, Ensembl id ENSG00000111666, with SNP rs10860794). Visual examination of the expression levels in each genotype class (Figure 5A), suggest that this SNP is an eQTL in all three tissues, with similar effects in each tissue (Figure 5B). This is supported by the joint analysis, which shows strong evidence for an eQTL , and assigns probability effectively 1 to the consistent configuration . However, as shown by the , at an FDR of 0.05, the tissue-by-tissue analysis calls the eQTL only in fibroblasts.

Given the disagreement between the results from our novel framework and the original analyses of these data, we checked the plausibility of our results by applying a previously-used method for examining pairs of tissues to these data ([13]). This analysis takes the best eQTL in each gene identified in one tissue, and then estimates the proportion of these () that are also eQTLs in a second tissue (by applying Storey's method [31] to their nominal in that second tissue, uncorrected for multiple comparisons). Unlike the tissue-by-tissue analysis above, this approach avoids thresholding of , and makes some allowance for incomplete power. However, unlike our framework, this approach can only be applied to compare pairs of tissues. Applying this approach to each pair yielded a mean estimate of (range to ), broadly consistent with our qualitative conclusion that a substantial proportion of eQTLs are shared among tissues.

Discussion

In this work, we have presented a statistical framework for analyzing and identifying eQTLs, combining data from multiple tissues. Our approach considers a range of alternative models, one for each possible configuration of eQTL sharing among tissues. We compute Bayes Factors that quantify the support in the data for each possible configuration, and these are used both to develop powerful test statistics for detecting genes that have an eQTL in at least one tissue (by Bayesian model averaging across configurations), and to identify the tissue(s) in which these eQTLs are active (by comparing the Bayes factors for different configurations against one another). Our framework allows for heterogeneity of eQTL effects among tissues in which the eQTL is active, for different variances of gene expression measurements in each tissue, and for intra-individual correlations that may exist due to samples being obtained from the same individuals. For eQTL detection, our framework provides consistent, and sometimes substantial, gains in power compared to a tissue-by-tissue analysis and ANOVA or simple linear regression. Concerning the tissue specificity of eQTLs, our framework efficiently borrows information across genes to estimate configuration proportions, and then uses these estimates to assess the evidence for each possible configuration. When re-analyzing the gene expression levels in three cell types from 75 individuals ([12]), we found that there appears to be a substantial amount of sharing of eQTLs among tissues, substantially more than suggested by the original analysis.

In the next few years, we expect that expression data will be available on large numbers of diverse tissue types in sufficient sample sizes to allow eQTLs to be mapped effectively (for example, the NIH GTEx project aims to collect such data). The methods presented here represent a substantive step towards improved analyses that fully exploit the richness of these kinds of data. However, we also see several directions for potential extensions and improvements. First, our current framework can only partially deal with the challenges of large numbers of tissues. Specifically, because with tissues, there are possible configurations of eQTL sharing among tissues, some of our current methods, which consider all possible configurations, will become impractical for moderate (speculatively, above about 10, perhaps). Our test statistic partially addresses this problem, by allowing for heterogeneity while averaging over only configurations, which is practical for very large . Our simulation results suggest that is a powerful test statistic for identifying SNPs that are an eQTL in at least one tissue. However our preferred approach for identifying which tissues such SNPs are active in involves a hierarchical model that estimates the frequency of different patterns of sharing from the data, and this hierarchical model scales poorly with . In particular, having a separate parameter for each possible configuration is unattractive (both statistically and computationally) for large , and alternative approaches will likely be required. There are several possible ways forward here: for example, one would be to reduce the number of distinct configurations by clustering “similar” configurations together; another would be to focus less on the discrete configurations, and instead to focus on modeling heterogeneity in effect sizes in a continuous way - perhaps using a mixtures of multivariate normal distributions with more complex covariance structures than we allow here. We expect this to remain an area of active research in the coming years, especially since these types of issues will likely arise in many genomics applications involving multiple cell types, and not only in eQTL mapping.

Another important issue to address is that most future expression data sets will likely be collected by RNA-seq, which provides count data that are not normally distributed. Previous eQTL analyses of RNA-seq (e.g. [32]) have nonetheless performed eQTL mapping using a normal model, by first transforming (normalized) count data at each gene to the quantiles of a standard normal distribution. Although this approach would not be attractive in experiments with small sample sizes, with the moderate to large sample sizes typically used in eQTL mapping experiments this approach works well. As a first step, this approach could also be used to apply our methods to count data. However, ultimately it would seem preferable to replace the normal model with a model that is better adapted to count-based data, perhaps a quasi-Poisson generalized linear model ([33]); Bayes Factors under these models could be approximated using Laplace approximations, similar to the approximations used here for the normal model [21]. The quasi-Poisson model has the advantage over the normal transformation approach that it preserves the fact that there is more information about eQTL effects in tissues where a gene is high expressed than in tissues where it is low expressed. This information is lost by normal transformation. In our primary analyses here we addressed this by analyzing only genes that were robustly expressed in all tissues, but this is sub-optimal, and will become increasingly unattractive as the number of tissues grows.

Our analyses here assess (cis) eQTL sharing among tissues by performing association testing at the level of individual SNPs. A different approach to investigating eQTL sharing among tissues is to study the “cross-heritability” of expression levels among tissues (e.g. [34], [35]). These methods are based on polygenic models, and attempt to estimate the combined influence of all shared eQTLs; this contrasts with our analysis, where the focus is on sharing of individually-identifiable eQTLs of moderate-to-large effect. Both [34] and [35] estimate cross-tissue heritability to be low. [34], studying expression in Blood and Adipose tissues from Icelanders, estimated cross-tissue heritability as 3%; [35] obtained an estimate of mean genetic correlation close to zero for Blood and LCLs in monozygotic twins (). These results may appear to conflict with our results (both from our model-based approach, and the less-model-based pairwise analysis approach from [13]), which suggest that most large-effect cis eQTLs are shared among fibroblasts, LCLs and T cells. However, these low estimates of cross-tissue heritability reflect not only the extent of sharing of eQTLs, but also the absolute size of the eQTL effects. If eQTL effects are small, explaining only a small proportion of the total variance in gene expression, then cross-tissue heritability will be also small, even if all eQTLs have exactly the same effect in all tissues. Thus, to assess eQTL sharing in the heritability-based approaches, it is helpful to contrast cross-tissue heritability, , with within-tissue heritability, , (which is also affected by eQTL effect size, but not by sharing). Specifically, within the polygenic model it can be shown that the correlation coefficient of the eQTL effects in two tissues and is: . Applying this to the cis estimates of and from [34], for adipose and blood, yields . Although this estimate of effect correlation within a polygenic model, is not directly comparable with our estimate of sharing of eQTLs in a decidedly non-polygenic model (and for different cell types!), this result suggests that the two analyses may be less in conflict than they initially appear.

Methods

Software implementing our methods are available on the website http://stephenslab.uchicago.edu/software.html.

Bayesian Methods for Mapping Multiple-Tissue eQTLs

Models for multiple-tissue eQTLs

For each tissue, we model the potential genetic association between a target SNP and the expression levels of a target gene by the simple linear regression model (1). In vector form, this model is represented by

where indexes one of the tissue types examined and the vectors and denote the expression levels, the genotypes of the samples and the residual errors respectively for the tissue type. The intercept term, , and the residual error variance, are allowed to vary with tissue type. The regression coefficient denotes the effect of the eQTL in tissue , but we follow [21], [28] in using the (unitless) standardized regression coefficient , as the main measure of effect size. As a result, inference is invariant to scale transformations of the response variables () within each tissue.

When the tissue samples are taken from the same individuals we allow that the observations on the same individual may be correlated with one another. Specifically, let denote the matrix of residual errors, we assume it to follow a matrix-variate normal (MN) distribution, i.e.,

That is, the vectors are independent and identically distributed as . The (unknown) covariance matrix quantifies the correlations between the tissues; it can vary from gene to gene and is estimated from the data (see below). [When the tissue samples are collected from different individuals then we assume their error terms are independent; methods for this case are given in [21].

Prior on effect sizes

A key component of our Bayesian model is the distribution , where denotes hyper-parameters that are to be specified or estimated from the data. (In the main text we used to simplify exposition, but we actually work with the standardized effects .) Of course, if then by definition. So it remains to specify the distribution of the remaining values for which .

We use the distribution from [21] (see also [19], [20]), which provides a flexible way to model the heterogeneity of genetic effects of an eQTL in multiple tissues. Specifically, [21] consider a distribution , with two hyper-parameters, , in which the non-zero effects are normally distributed about some mean , which itself is normally distributed:

and

Note that controls the variance (and hence the expected absolute size) of , and controls the heterogeneity (indeed, is the correlation of for different subgroups ). If then this model corresponds to the “fixed effects” model in which the effects in all subgroups are equal (e.g. [23]).

To allow for different levels of effect size and heterogeneity, [21] use a fixed grid of values , with the th grid point having weight . Thus

In all our applications here we consider the grid of values fixed, and treat the weights as hyper-parameters (so ), which can be either fixed or estimated.

Choice of grid for

We define a grid of values for by specifying a set of values for the average effect size, , and a set of values for the heterogeneity , and then taking the grid to be all possible combinations of values. For all methods we use , which is designed to span a wide range of eQTL effect sizes (see Text S1 as well as Figure S2). For and we allow for only a limited range of heterogeneity: . In this way we assume that when the eQTL is present in multiple tissues, it has a similar (but not necessarily identical) effect in each tissue. For we allow a much wider range of heterogeneity: . The rationale here is that the large heterogeneity values will help capture eQTLs that are present in only a subset of tissues, a feature that is not otherwise captured by as it averages over a small number of configurations.

Choice of weights and

Let denote the number of elements of that are equal to 1 (i.e. the number of tissues in which the eQTL is active in configuration ), and recall that .

For we fix the weights so that they put weight on all possible non-zero values for (), and, conditional on , put equal weight on all configurations with that value for . Thus . In addition we fix the grid weights to be equal on all grid values.

For we put non-zero weights on only the consistent configuration () and configurations with an eQTL in a single tissue (). We set so that it puts weight 0.5 on each of and . Conditional on we assume all possibilities are equally likely. Thus if and 0.5/S if . Again, we fix the grid weights to be equal on all grid values (but with the larger grid for heterogeneity described above).

For we estimate the weights, from the data using a hierarchical model to combine information across genes, as described below.

Bayes Factor computation

To complete model specification, we use (limiting, diffuse) prior distributions for the nuisance parameters and , as in [36]. Under these priors we can compute the Bayes Factor in (2) using

where is the total number of grid points and is given by

where and denote the collection of expression levels and genotypes for a target gene-SNP pair across all tissue types respectively. We use analytic approximations for these Bayes Factors based on Laplace approximation, given in [21], . In particular, we use the approximation which in functional forms is connected to Frequentist's score statistic.

Bayesian Hierarchical Model

For we use a hierarchical model, similar to [24], [25], which combines information across genes, to estimate the grid weights 's and configuration weights 's. Following both [25], [26] we make the simplifying assumption that each gene has at most one eQTL (which may be active in multiple tissues), and that each SNP is equally likely to be the eQTL. Let be the number of SNPs in the cis-region for gene . Then, if denotes the Bayes Factor (12) computed for SNP in gene , the “overall Bayes Factor” measuring the evidence for an eQTL in gene , , is obtained by averaging over the possible eQTL SNPs, the possible configurations , and the grid of values for , weighting by their probabilities:

Furthermore, if we let denote the probability that each gene follows the null (i.e. contains no eQTL) then the likelihood for gene , as a function of , is given by

The overall likelihood for our hierarchical model is obtained by multiplying these likelihoods across genes:

Note that although the expression levels for different genes are not independent, because the SNPs being tested in different genes are mostly independent this independence assumption for the likelihoods across genes is a reasonable starting point. We have developed an EM algorithm to estimate the parameters by maximum likelihood (see Supplementary information).

Relaxation of “one cis-eQTL per gene” assumption

To relax the “one cis-eQTL per gene” assumption we adopt the following procedure. First we compute the posterior probability of each SNP being the sole eQTL for each gene (i.e. only allowing one cis-eQTL per gene) with a set of default parameters, and use these to identify the top SNP for each gene (i.e. the one with the largest posterior probability of being the eQTL).

For each gene, separately in each tissue, we compute the residuals of its expression level after regressing out the effect of the top SNP. If these residuals are strongly associated with a SNP then this is evidence for that SNP being a second independent eQTL for that gene. Therefore, to allow for potentially more than one eQTL per gene we treat these residuals as defining a second set of “artificial” expression data for each gene and each tissue, and fit the hierarchical model using both the original and the artificial expression data.

Simulation Procedures

For our simulations, when simulating SNP-gene pairs, the genotypes at each SNP in each individual were simulated as Binomial(2,0.3): that is, with minor allele frequency 30% and assuming Hardy-Weinberg equilibrium. Phenotypes with eQTLs were simulated, with effect size based on an expected proportion of variance explained (PVE) of 20%; (see Text S1). For Figure 1A and 1B the error variances (one per tissue) were all equal to 1. For Figure 1C the error variances were randomly drawn from , all equally likely.

The ANOVA/LR Method

The ANOVA/LR method uses the same linear model as our Bayesian methods (1), except that the residual errors are assumed to be equal across tissues . Within this model we tested the global null hypothesis ( for all ) using an test comparing the null model with the unconstrained alternative ( unconstrained). See Text S1.

Preprocessing of the Dataset from Dimas et al.

The phenotypes from Dimas et al. ([12]) were retrieved from the Gene Expression Omnibus (GSE17080). We mapped the 22,651 non-redundant probes to the hg19 human genome reference sequence (only the autosomes) using BWA ([37]), kept 19,965 probes mapping uniquely with at most one mismatch, and removed the probes overlapping several genes from Ensembl. This gave us 12,046 genes overlapped by 16,453 probes. For genes overlapped by multiple probes, we chose a single probe at random. In our analyses we considered only genes that were robustly expressed in all tissues. A gene was considered robustly expressed in a given tissue if its mean expression level across individuals in this tissue was larger than or equal to the median expression level of all genes across all individuals in this tissue. As a result, we focused on 5012 genes.

Genotypes were obtained from the European Genome-phenome Archive (EGAD00000000027). We extracted the genotypes corresponding to the 85 individuals for which we had phenotypes and converted the SNP coordinates to the hg19 reference using liftOver ([38]). To detect outliers, we performed a PCA of these genotypes using individuals from the CEU, CHB, JPT and YRI populations of the HapMap project using EIGENSOFT ([39]). As in the original study, we identified 10 outliers and removed them from all further analyses, which were therefore performed on 75 individuals.

Gene expression measurements suffer from various confounders, many of which may be unmeasured ([40]), but which can be corrected for using methods such as principal components analysis (PCA). Following [32], we applied PCA in each tissue separately on the 501275 matrix of expression levels of each gene in each individual. We sorted principal components (PCs) according to the proportion of variation in the original matrix they explain, and selected PCs so that adding another PC would explain less than 0.0025% of the variation. As a result, this procedure identified 16 PCs in Fibroblasts, 7 in LCLs and 15 in T-cells. We then regressed out these PCs from the original matrix of gene expression levels, and used the residuals as phenotypes for all analyses.

All methods we compared assume that the errors are distributed according to a Normal distribution. Before analysis we therefore rank-transformed the expression levels at each gene to the quantiles of a standard Normal distribution ([28]).

Permutation Procedures

On the data set from Dimas et al., we assessed the performance of two methods, the tissue-by-tissue analysis and the BMA joint analysis, by comparing the number of genes identified as having at least one eQTL in any tissue, at a given FDR. For each method, we defined a test statistic, which was computed for each gene. For the tissue-by-tissue analysis, the test statistic is the minimum of the linear regressions between the given gene and each cis SNP in each tissue (so the minimum is taken across all SNPs and all tissues). For the BMA joint analysis, the test statistic is the average of the Bayes Factors for the given gene and each cis SNP. (When applying the tissue-by-tissue analysis to test for eQTLs in a single tissue, the test statistic is the minimum of the linear regressions between the given gene and each cis SNP in that tissue.)

In each case we converted the test statistic to a for each gene, testing the null hypothesis that the gene contains no eQTL in any tissue, by comparing the observed test statistic with the value of the test statistic obtained on permuted data obtained by permuting the individuals labels (using the same permutations in each tissue to preserve any intra-individual correlations between gene expression in different tissues). Specifically, let denote the total number of permutations (we used ), the value of the test statistic for gene on the non-permuted data, and the value of the test statistic on the -permuted data. The for gene from the tissue-by-tissue analysis is: . For the BMA joint analysis, the is: . Note that permutations were performed for each gene, since the null distribution of the test statistic will vary across genes (not least because the genes have different numbers of SNPs in their cis candidate region; see Figure S3).

From the calculated for each gene we estimate using the qvalue package ([31]), and determine the number of genes having at least one eQTL in any tissue at an FDR of by computing the number of genes with .

When performing the tissue-by-tissue analysis on a single tissue, we performed the permutations in each tissue separately.

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4