PaperlayerView exact revision on arXiv

Source-linked AI summary

A statistical framework for joint eQTL analysis in multiple tissues

Timothée Flutre, Xiaoquan Wen, Jonathan Pritchard, Matthew Stephens

arXiv:1212.4786v1q-bio.QMq-bio.GNstat.AP

TL;DR

The paper addresses the challenge of analyzing eQTLs across multiple tissues when effects may be shared in some tissues and inactive in others. It develops a configuration-based framework that jointly models these patterns and applies it to three tissues. The framework identifies 63% more eQTLs than tissue-by-tissue analysis at FDR=0.05, while the results suggest that most detectable eQTLs are shared across all three tissues.

  • Problem

    Existing tissue-by-tissue analyses do not fully exploit shared eQTLs across tissues or formally account for incomplete power when comparing overlaps.

  • Method

    The framework evaluates alternative binary configurations of eQTL activity across tissues, using Bayesian model averaging and Bayes Factors to detect eQTLs and identify active tissues.

  • Results

    63% more eQTLs were identified by BFBMA than by tissue-by-tissue analysis at FDR=0.05, with 1022 versus 627 eQTLs.

  • Takeaways & Limitations

    The framework increases power for joint eQTL mapping and supports estimating tissue-sharing patterns while allowing inactive tissues and heterogeneous effects.

  • Takeaways & Limitations

    Considering all possible sharing configurations becomes challenging as the number of tissues grows because S tissues yield 2^S configurations.

Abstract

from arXiv · show

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

Understanding regulatory variation across tissues is important for interpreting phenotypes, but existing analyses do not fully exploit shared and tissue-specific eQTL patterns. The paper introduces a joint framework that combines information across tissues while allowing heterogeneity and estimates patterns of sharing.

  • Motivation: Regulatory variation across tissues may help connect genetic variants to disease biology and illuminate tissue differentiation.The motivation includes identifying which tissues are affected by regulatory variants associated with GWAS hits.
  • Limitations of existing analyses: Tissue-by-tissue analyses fail to use commonalities among tissues to improve power for detecting shared eQTLs.Existing tools are also limited in jointly analyzing all tissues while allowing differences among eQTLs.
  • Related approaches: Joint analyses can increase power for eQTLs with similar effects across tissues, but existing ANOVA and weighted Z-score approaches have different advantages.The cited methods provide greater power than tissue-by-tissue analysis, while differing in how they address effect heterogeneity.
  • Proposed framework: The proposed framework combines multiple-tissue information in a hierarchical model that estimates the relative frequency of eQTL-sharing patterns across genes.It integrates GWAS meta-analysis methods that allow heterogeneous effects among groups.
  • Proposed framework: The framework models heterogeneity across several tissues, allows tissue-specific variances and intra-individual correlations, and explicitly estimates which tissues share eQTLs.It can also be applied across experimental platforms, datasets, or populations rather than only tissue types.
  • Study design: The study uses simulations and data from fibroblasts, LCLs, and T-cells to evaluate power and tissue-consistent eQTLs.The authors report a large power gain over tissue-by-tissue analysis and a higher rate of tissue-consistent eQTLs than previous analyses.

Results

The framework models tissue-specific activity, heterogeneity, and sharing patterns while combining evidence across genes and tissues. Simulations and applications show increased power over tissue-by-tissue analysis, especially for shared eQTLs, and support substantial cross-tissue sharing.

  • Framework: A configuration encodes whether an eQTL is active or inactive in each tissue, while P(β|γ, θ) allows effect sizes and cross-tissue heterogeneity to vary.The global null is the configuration in which the SNP is inactive in every tissue.
  • Simulation results: Joint analyses outperform tissue-by-tissue analysis for eQTLs occurring in multiple tissues, with larger gains as the number of sharing tissues increases.Tissue-by-tissue analysis performs best for single-tissue eQTLs, but only slightly better than joint approaches in that setting.
  • Framework: Bayesian Model Averaging averages over possible eQTL configurations, combining advantages of tissue-by-tissue and joint analyses across diverse eQTL patterns.This strategy is designed to detect eQTLs present in either single or multiple tissues.
  • Application results: Joint analyses outperformed tissue-by-tissue analysis for every tissue, while BFHM BMA outperformed BFBMA by learning sharing patterns from the data.The gain was greater for Tissue 1 than Tissue 2, illustrating larger benefits for tissues with smaller sample sizes.
  • Application results: BFBMA identified 1022 eQTLs at FDR=0.05, 63% more than the 627 identified by tissue-by-tissue analysis at the same FDR.BFBMA also detected 94% of the eQTLs identified by tissue-by-tissue analysis.
  • Sharing estimates: The hierarchical model estimated 8% of eQTLs as specific to one tissue and 88% as common to all three tissues, with a 95% CI of 84%- 93%.Pairwise π1 estimates averaged approximately 88% across tissue pairs, broadly supporting substantial sharing.

Discussion

The framework advances joint eQTL analysis while highlighting scalability and RNA-seq modeling challenges that motivate future extensions.

  • Discussion: The framework models eQTL-sharing configurations across tissues and uses Bayes Factors to support detection and model averaging.It considers one model for each possible sharing configuration and uses the resulting Bayes Factors to construct detection statistics.
  • Discussion: For S tissues, evaluating all 2^S sharing configurations can become impractical above about 10 tissues.The BFBMAlite statistic addresses this by averaging over S + 1 configurations while allowing heterogeneity.
  • Discussion: These challenges are expected to recur in genomics applications involving multiple cell types beyond eQTL mapping.The authors identify this as an area for continued research as datasets expand across diverse tissues.
  • Discussion: RNA-seq count data violate the normality assumption used by the current framework, although transformed counts can provide an interim application for moderate to large samples.A future quasi-Poisson generalized linear model is suggested as a better-adapted alternative.

Materials and Methods

The framework models tissue-specific eQTL associations while allowing active effects to be shared, heterogeneous, or absent across tissues. Bayesian effect-size and configuration models combine evidence across tissues and genes, with tissue correlations accommodated for shared samples.

  • Tissue-level model: The tissue-level model regresses expression on genotype, allowing tissue-specific intercepts, residual variances, and eQTL effects.Standardized effects b_s are used so inference is invariant to within-tissue response scaling.
  • Tissue-level model: Correlated residual errors across tissues are modeled when samples come from the same individuals, using a gene-specific covariance matrix.Samples from different individuals are instead treated as having independent error terms.
  • Prior on effect sizes: The prior sets inactive tissue effects to zero and uses a flexible distribution for active effects, controlling typical effect size and cross-tissue heterogeneity.A hierarchical normal construction models nonzero effects around a shared mean, while variance components determine heterogeneity and correlation.
  • Prior on effect sizes: The effect-size grid spans five average-effect values and uses narrower heterogeneity settings for BFBMA and BFHM than for BFBMAlite.BFBMA and BFHM use H = {0, 0.25}, whereas BFBMAlite uses H = {0, 0.25, 0.5, 0.75, 1}.
  • Sharing configurations: Configuration weights encode alternative patterns of eQTL activity, including equal weighting across the numbers of active tissues and a lite model emphasizing single-tissue and fully shared configurations.For BFBMAlite, each of |γ| = 1 and |γ| = S receives weight 0.5.
  • Bayes factors and hierarchical modeling: Bayes factors integrate expression, genotype, nuisance parameters, effect-size configurations, and grid values, using Laplace approximations connected functionally to a frequentist score statistic.The hierarchical model averages evidence over candidate SNPs, tissue configurations, and effect-size grid points.
  • Bayes factors and hierarchical modeling: A hierarchical likelihood combines information across genes to estimate null, configuration, and grid weights, while initially assuming at most one cis-eQTL per gene.The model treats this independence across genes as a reasonable starting point because SNPs tested in different genes are mostly independent.
  • Model extension: To relax the one-eQTL assumption, the procedure first estimates posterior probabilities under that restriction and then identifies each gene’s top SNP.The top SNP is the one with the largest posterior probability of being the eQTL.

A Computational algorithm for fitting hierarchical model

The hierarchical model estimates its parameter set by maximum likelihood using an expectation-maximization algorithm.

  • Algorithm: The EM algorithm infers Θ = (π0, η, λ) by maximum likelihood for the hierarchical model.The parameter set includes the null probability and model-weight parameters.

A.1 Notations

The computational formulation represents gene-level, SNP-level, tissue-configuration, and effect-size uncertainty with latent variables. EM alternates posterior expectation and parameter maximization, with convergence monitoring and profile-likelihood intervals.

  • Latent variables: A latent indicator z_k records whether gene k has any eQTL in its cis-region.
  • Latent variables: The one-cis-eQTL assumption represents the true eQTL SNP with s_k, which has at most one active entry.
  • Latent variables: Configuration indicators c_kp encode which of the 2^S − 1 tissue-activity patterns applies to each gene–SNP pair.
  • Latent variables: Effect-size indicators w_kp select the prior effect-size grid for active tissues, and the resulting indicators form W_k.
  • EM fitting: The EM algorithm treats z_k, s_k, c_k, and w_k as missing data while iterating expectation and maximization steps.Initialization uses parameter values and iterations stop when successive log-likelihood increases become sufficiently small.
  • Uncertainty assessment: Profile-likelihood confidence sets are constructed for estimated parameters such as π0.

B.1 Simulate eQTL data via the proportion of variance explained

The simulation framework generates multi-tissue eQTL phenotypes by specifying variance explained, allele frequency, and effect heterogeneity. Simulated data can then be tested with an ANOVA/likelihood-ratio model.

  • Simulation parameterization: For each gene–SNP pair, simulations generate data across S tissues with effects active in selected tissues.The phenotype variance explained by genotype is used to parameterize the simulation.
  • Simulation parameterization: The approximation h = (φ2 + ω2) × 2f(1 − f) / ((φ2 + ω2) × 2f(1 − f) + 1) links effect-size variance and minor allele frequency to explained variance.
  • Simulation parameterization: Fixing h and minor allele frequency determines the effect-size variance, while fixing heterogeneity separates φ2 and ω2.The resulting parameters generate the shared mean effect and tissue-specific effects.
  • Phenotype generation: Expression values are simulated for individuals and tissues, with residual scale optionally fixed at σs = 1.
  • Testing simulated data: All tissue expression values and repeated genotypes are assembled into vectors for an ANOVA/likelihood-ratio test of genotype effects.The implementation compares models with tissue indicators alone versus tissue-by-genotype interactions.

B.3 Calculate the empirical FDR from simulated eQTL data

The empirical FDR framework classifies simulated gene-SNP pairs by discovery status and known truth, then fixes the empirical FDR at 5% to calculate corresponding TPRs. This enables comparison of methods through their TPR and FPR values.

  • For the tissue-by-tissue method, the test statistic is the minimum P-value across tissues; Bayesian analysis uses the Bayes Factor, while ANCOVA uses the genotype-effect interaction P-value.
  • Simulated gene-SNP pairs are classified by whether they are called eQTLs and whether they are true eQTLs.
  • The empirical FDR is fixed at 5% to determine a test-statistic cutoff and calculate the resulting true positive rate.Because simulations identify true eQTLs, true positives and false positives can be counted directly.
  • An iterative procedure sorts test statistics or Bayes Factors, evaluates called eQTLs, and stops when the empirical FDR reaches the 5% threshold.
  • At a common empirical FDR of 5%, methods can be compared because their true positive and false positive rates differ.

C.1 Hierarchical model fed with Bayes Factors from residuals

The hierarchical model is initialized with Bayes Factors computed for each gene-SNP pair, then rerun on residual phenotypes after regressing out each gene’s best SNP. Under the at-most-one-eQTL assumption, expected parameter patterns indicate weak effects and substantial uncertainty.

  • Bayes Factors are computed across grid values and configurations for each gene-SNP pair, followed by residualization against each gene’s best SNP.
  • The hierarchical model is then launched using Bayes Factors computed from the residual phenotypes.
  • If each gene has at most one eQTL, the estimated π0 should be very high, the smallest grid value most probable, and configuration credible intervals very large.These patterns correspond to most genes having no eQTL, very small effect sizes, and high uncertainty.

C.2 Configuration proportions from all genes without removing expression PCs

The analysis was repeated across all 12,046 genes without selecting genes robustly expressed in all tissues or removing expression PCs. The resulting configuration proportions were qualitatively similar to the main analysis.

  • The analysis included all 12,046 genes without pre-selecting genes robustly expressed in all three tissues.
  • Expression PCs were not removed in this all-gene analysis.
  • The configuration proportions estimated by the EM algorithm were qualitatively similar to those from the filtered, PC-adjusted analysis.
Loading 1212.4786v1…