New insights into the genetic control of gene expression using a Bayesian multi-tissue approach

Abstract

The majority of expression quantitative trait locus (eQTL) studies have been carried out in single tissues or cell types, using methods that ignore information shared across tissues. Although global analysis of RNA expression in multiple tissues is now feasible, few integrated statistical frameworks for joint analysis of gene expression across tissues combined with simultaneous analysis of multiple genetic variants have been developed to date. Here, we propose Sparse Bayesian Regression models for mapping eQTLs within individual tissues and simultaneously across tissues. Testing these on a set of 2,000 genes in four tissues, we demonstrate that our methods are more powerful than traditional approaches in revealing the true complexity of the eQTL landscape at the systems-level. Highlighting the power of our method, we identified a two-eQTL model (cis/trans) for the Hopx gene that was experimentally validated and was not detected by conventional approaches. We showed common genetic regulation of gene expression across four tissues for approximately 27% of transcripts, providing >5 fold increase in eQTLs detection when compared with single tissue analyses at 5% FDR level. These findings provide a new opportunity to uncover complex genetic regulatory mechanisms controlling global gene expression while the generality of our modelling approach makes it adaptable to other model systems and humans, with broad application to analysis of multiple intermediate and whole-body phenotypes.

Figure 1. Evidence of pleiotropic eQTLs detected by the SBMR model.

(A, E) For each gene, the set of markers associated with high marginal posterior probability of inclusion corresponds to the filtered best model found by SBMR, showing monogenic control for Cd36 gene (marker Cd36) and polygenic control for Ascl3 gene (markers D1Rat55 and D7Mit8). The marginal probability of inclusion is calculated conditionally on all visited models whose Bayes Factor is above the calibrated threshold at 5% FDR level (see Materials and Methods). (B, C, F, G, H) Systemic effects of pleiotropic eQTLs detected by the SBMR model. For each gene we report the raw empirical correlation of gene expression across four tissues and the posterior mean of the residual correlation matrix. Posterior correlations in panels C, H were simulated conditionally on the filtered best model that coincide with the noticeable effects (see Materials and Methods): marker Cd36 for Cd36 gene and markers D1Rat55 and D7Mit8 for Ascl3 gene, respectively. In panel G, the posterior correlations were generated conditionally on the cis-eQTL only (marker D1Rat55). The size of the correlation is colour coded and reported in each graph. (D, I, L) Box-plots of the posterior density of the effect size (see Text S1) for the eQTLs with noticeable effect are reported for each tissue. These illustrate the tissue-specific contribution provided by each eQTL to the pleiotropic effect. Tissues: A, adrenal; F, fat; H, heart; K, kidney.

Figure 2. Sensitivity and specificity of SBMR and alternative approaches.

Log-scale Receiver Operating Characteristic (ROC) curves of SMBR (blue), Hotelling's -test (red) and GFlasso (green) methods, using simulated data generated under five different scenarios. The scenarios for the pleiotropic eQTL are as follows: (A) one cis-eQTL; (B) one cis- and one trans-eQTLs; (C) two trans-eQTLs; (D) one cis-eQTLs and four trans-eQTLs and (E) four trans-eQTLs. In each case we simulated strong (left panel), medium (central panel) and weak (right panel) correlation pattern among gene expression traits (see Text S1 for details).

Figure 3. Validation of polygenic regulation for Hopx.

The filtered best model for the regulation of Hopx indicates polygenic control of gene expression by two co-existing eQTLs, D14Rat36 and D2Rat136. The marginal posterior probability of inclusion for the cis- (D14Rat36) and trans-eQTL (D2Rat136) is reported in panel (A). RT-PCR data showing relative Hopx expression in the BXH/HXB RI strains by BN and SHR genotypes at peak markers D14Rat36 (left panel) and D2Rat136 (right panel), (B). The cis-eQTL is identified by all methods (SSM: FDR <5%; QTL Reaper: genome-wide corrected p-value, , FDR <5%), while the weaker trans-eQTL at marker D2Rat136 (indicated by an arrow) is not significantly detected by either the SSM (panel C) or QTL Reaper (panel D) methods. This shows the power of the SBR model to identify both small (trans-acting) and big (cis-acting) genetic effects that can simultaneously determine variation in gene expression. Relative expressions are reported as mean ± sem. (, ).