Source-linked AI summary
Integrative modeling of eQTLs and cis-regulatory elements suggest mechanisms underlying cell type specificity of eQTLs
Christopher D Brown, Lara M Mangravite, Barbara E Engelhardt
TL;DR
The paper addresses the limited understanding of eQTL genetic architecture and cell type specificity needed to interpret GWAS associations mechanistically. It analyzes eQTLs across eleven studies and seven cell types, integrates them with ENCODE cis-regulatory-element data, and builds a classifier. The analyses find pervasive allelic heterogeneity, frequent cell type specificity, regulatory-element associations with replication, and successful prediction of eQTL specificity without target-cell expression data.
Problem
Interpreting GWAS associations with eQTLs requires better understanding of eQTL genetic architecture and cell type specificity.
Method
The study analyzes eleven eQTL studies across seven cell types, integrates eQTLs with ENCODE CRE data, and trains a random forest classifier using 526 CRE data sets.
Results
The analyses show pervasive cis-eQTL allelic heterogeneity, frequent cell type specificity, CRE-associated replication patterns, and successful prediction of specificity without target-cell expression data.
Takeaways & Limitations
Integrating eQTLs with cell type-specific CREs can support mechanistic interpretation and estimate eQTL replication between cell types when target-cell expression data are unavailable.
Abstract
from arXiv · showhide
Genetic variants in cis-regulatory elements or trans-acting regulators commonly influence the quantity and spatiotemporal distribution of gene transcription. Recent interest in expression quantitative trait locus (eQTL) mapping has paralleled the adoption of genome-wide association studies (GWAS) for the analysis of complex traits and disease in humans. Under the hypothesis that many GWAS associations tag non-coding SNPs with small effects, and that these SNPs exert phenotypic control by modifying gene expression, it has become common to interpret GWAS associations using eQTL data. To exploit the mechanistic interpretability of eQTL-GWAS comparisons, an improved understanding of the genetic architecture and cell type specificity of eQTLs is required. We address this need by performing an eQTL analysis in four parts: first we identified eQTLs from eleven studies on seven cell types; next we quantified cell type specific eQTLs across the studies; then we integrated eQTL data with cis-regulatory element (CRE) data sets from the ENCODE project; finally we built a classifier to predict cell type specific eQTLs. Consistent with prior studies, we demonstrate that allelic heterogeneity is pervasive at cis-eQTLs and that cis-eQTLs are often cell type specific. Within and between cell type eQTL replication is associated with eQTL SNP overlap with hundreds of cell type specific CRE element classes, including enhancer, promoter, and repressive chromatin marks, regions of open chromatin, and many classes of DNA binding proteins. Using a random forest classifier including 526 CRE data sets as features, we successfully predict the cell type specificity of eQTL SNPs in the absence of gene expression data from the cell type of interest. We anticipate that such integrative, predictive modeling will improve our ability to understand the mechanistic basis of human complex phenotypic variation.
Introduction
The study addresses the difficulty of interpreting eQTLs across cell types by quantifying their specificity and integrating eQTLs with regulatory-element data. It analyzes eQTLs across multiple studies and cell types, then develops a classifier to predict cell type specificity without expression data from the target cell type.
- Motivation: The cell type specificity of eQTLs complicates mechanistic interpretation of GWAS associations when the phenotype-relevant cell type is unknown or poorly represented in available eQTL data.An eQTL observed in lymphoblastoid cell lines may not produce the same molecular phenotype in the relevant tissue.
- Approach: The study analyzes eQTL data from eleven studies spanning seven cell types to quantify cell type-specific eQTL patterns.
- Approach: It integrates eQTL SNPs with 526 cis-regulatory-element data sets, including cell type-specific transcription-factor binding and open-chromatin data.
- Approach: The study builds a random forest classifier to predict eQTL cell type specificity without additional gene expression data from the cell type of interest.
Results
Across eleven studies spanning seven cell types, cis-eQTLs showed pervasive allelic heterogeneity, cell type-specific replication, and spatially structured overlap with regulatory elements. Replication and regulatory-element overlap varied with study power, SNP location, eQTL tier, and cell type.
- Allelic heterogeneity: 29% of eQTL-regulated genes were independently associated with at least two SNPs in at least one study.Within studies, the fraction ranged from 3–22% at FDR ≤5%.
- Cell type specificity: Approximately 80% of eQTLs replicated within the same cell type, compared with approximately 60% between different cell types.Replication increased with discovery significance and decreased with distance from the TSS and eQTL tier.
- Cell type specificity: Distal eQTLs were less reproducible than promoter-proximal eQTLs, including across studies of the same cell type.This location effect must be considered when comparing between-cell-type replication.
- Cell type specificity: Primary eQTL SNPs were more reproducible and less cell type specific than additional independently associated SNPs.For CAP LCL, 63.4% of primary and 73.0% of secondary SNPs were cell type specific.
- CRE overlap: eQTL SNPs were enriched near active regulatory elements and depleted in repressed chromatin, heterochromatin, and regions with intervening CTCF sites.DHS overlap was greatest adjacent to the TSS, whereas negative regulatory overlap was typically lowest from immediately upstream of the TSS through the gene body.
- CRE overlap: Primary and secondary CAP LCL eQTL SNPs were associated with 134/166 and 100/166 LCL CRE classes, respectively.Independently associated SNPs less than 20 kb apart were more than twice as likely as background cis-SNPs to have an intervening insulator: 55.7% versus 22.6%.
- CRE overlap: Reproducible eQTLs were more likely to overlap active CRE classes and less likely to overlap repressed CRE classes than non-reproducible eQTLs.Associations with reproducibility were observed across 20/164 LCL CRE data sets and 31/150 HepG2 CRE data sets.
Discussion
The analyses show that cis-eQTL replication differs by cell type and regulatory context, while revealing methodological limits and practical uses for interpreting GWAS associations.
- Replication patterns: Within-cell-type replication is consistently higher than between-cell-type replication, indicating that many cis-eQTLs are cell type specific.Non-primary eQTLs also replicate substantially and are more often promoter distal and cell type specific.
- Prediction and application: A random forest classifier successfully predicted eQTL cell type specificity from independently derived CRE data, although substantial improvement remains possible.The predictions use eQTL location, cell types, known or predicted CREs, and gene function.
- Regulatory context: eQTLs overlap active CREs more and repressive CREs less when measured in the same cell type than across different cell types.Activating CRE overlap predicts independent replication more strongly when CREs and expression data come from the same cell type.
- Regulatory context: Cell type-specific CRE overlap is enriched among cell type-specific eQTLs, supporting distinct regulatory mechanisms for these associations.The analyses included enhancer, promoter, repressive chromatin, open-chromatin, and DNA-binding-protein CRE classes.
- Limitations: The analysis created multiple expression values per gene by clustering probes across platforms and combined cluster-level eQTL results only after analysis.Future work will identify eQTLs jointly across probe clusters.
- Limitations: CRE colocalization does not establish whether a CRE permits SNP functionality or mediates the regulatory effect, because causal directionality was not examined.The authors identify directionality as an area for current research.
Materials and Methods
The study assembled genotype data from multiple sources and applied standardized quality-control procedures before downstream analysis.
- Genotype quality control: Genotype data came from public databases or individual investigators and were filtered using call-rate and Hardy–Weinberg-equilibrium criteria.Individuals and SNPs with call rates below 90% were removed or classified as missing, and SNPs deviating from HWE at p < 1 × 10^-4 were removed.
- Study-specific processing: Merck liver genotypes were restricted to SNPs that were more than 90% unimputed, with imputed genotypes removed before applying identical filters.The procedure was intended to represent the original genotyping data.
- Study-specific processing: Harvard study data retained 540 cerebellum, 678 prefrontal-cortex, and 463 visual-cortex individuals after ungenotyped individuals were removed.The genotypes were matched to indexed individuals and filtered using the same criteria.
Genotype imputation
Genotype imputation expanded variant coverage, while probe alignment, expression filtering, normalization, and probe clustering prepared gene-expression traits for analysis.
- Genotype imputation: Genotypes were imputed with BIMBAM to the HapMap phase 2 CEPH 3.8 × 10^6 SNP set after removing variants with MAF below 0.01 or missing data.Caucasian-only studies used 60 unrelated CEPH references; the UChicago liver study also used 60 unrelated YRI references.
- Probe processing: Probes were retained when they uniquely aligned to the human reference genome or adequately aligned to a RefSeq transcript at at least 90% identity.Other probes were removed.
- Expression processing: Genes below study-appropriate expression thresholds were removed using negative controls, spike-ins, or observed mean–variance relationships.The thresholds defined when a gene was considered expressed.
- Expression processing: Expression data were cleaned, background-corrected, log2 transformed, and missing values imputed with k-nearest neighbors using k = 10.Flagged or poorly extracted features were treated as missing, and negative adjusted intensities were set to half the minimum positive array value.
- Probe clustering: Multiple probes per gene were clustered into at most min(4, p) groups, whose PC-corrected means served as separate expression proxies.Each probe cluster was modeled independently under the assumption that uncorrelated probe sets may represent independent transcripts or regulatory units.
eQTL mapping
The eQTL mapping model tested SNP–expression associations separately across genes, samples, and studies using Bayesian regression.
- Association mapping: Bayesian regression quantified associations between each SNP and residual gene-expression trait for each gene across samples within each study.The model averaged over plausible effect sizes for additive and dominant inheritance models and used mean-imputed genotypes except for Stranger LCLs.
Summarizing eQTLs
The analysis identifies eQTL associations while controlling false discovery through study-specific permutation tests. Results reported as eQTLs satisfy an FDR threshold of 5%.
- False-discovery control: FDR was estimated by permuting sample indices identically across genes and comparing original with permuted association counts at each log10 BF cutoff.A single complete permutation was performed for each study because of computational resource requirements.
- False-discovery control: Unless otherwise noted, reported eQTL associations were significant at FDR ≤5%.
Multivariate analysis
The study uses a two-step LD-based Bayesian approach to identify independently associated SNPs and compares it with forward stepwise regression. The comparison provides a more complete but still conservative estimate of allelic heterogeneity for a selected gene subset.
- LD-based Bayesian analysis: Approximately 20,000 conditional QTL scans were handled using a two-step approach that first selected highly associated SNPs by LD block, then recomputed conditional Bayes factors.Candidate SNPs were selected within 1 Mb of each gene’s transcription start or end site.
- Forward stepwise comparison: Forward stepwise regression added SNPs that most improved BIC until no further improvement occurred, using all SNPs within 1 Mb of gene boundaries.BIC accounts for model likelihood and parameter count to reduce overfitting.
- Forward stepwise comparison: The forward-selection comparison was performed on CAP LCL genes with significant allelic heterogeneity at FDR ≤5%.The approach is tractable but not exhaustive, yielding a more complete yet conservative allelic-heterogeneity estimate.
Analysis of Gene Ontology annotation enrichment
Gene Ontology and pathway enrichment analyses tested whether allelic-heterogeneity and replication categories were enriched for biological processes, molecular functions, or KEGG pathways.
- Enrichment analysis: DAVID enrichment analyses evaluated GO Biological Process, Molecular Function, and KEGG pathway terms using an FDR ≤20% threshold.The analyses examined CAP LCL, UChicago liver, and Harvard cerebellum discovery data sets.
- Enrichment analysis: Enrichment was assessed among allelic-heterogeneity genes and among genes that did or did not replicate within and between cell types.
Replication quantification
Replication was defined using evidence for the same SNP–gene pair in a target study, conditional on its observation in the replication cohort. This definition accounts for both Bayesian evidence and study-specific data availability.
- Replication definition: A SNP–gene association was considered replicated when the target-study log10 BF was at least 1.0 and the discovery association met its tier-specific FDR ≤5% cutoff.
- Replication definition: Replication frequency included only SNP–gene pairs observed in the relevant replication cohort.Pairs were excluded when the SNP failed quality control or the gene was absent from the microarray platform.
Comparison between eQTLs and functional genomic data sets
The study compares eQTL SNPs with cis-regulatory and genomic-feature data using logistic-regression models that control for positional, expression, association, and SNP-pair factors. It also tests whether regulatory-element overlap relates to replication within and between cell types.
- Overlap definition: SNP overlap included containment within a genomic element or location within 500bp of its boundary.This accounts for eQTL SNPs potentially tagging causal variants in linkage disequilibrium.
- CRE enrichment: Logistic regression modeled cis-regulatory-element overlap while controlling for SNP distance from the transcription start site and associated-gene expression level.The indicator variable distinguished observed eQTL SNPs from background cis-linked SNPs.
- Replication analysis: A second logistic-regression model estimated how cis-regulatory-element overlap relates to within-cell-type eQTL replication.The model controlled for SNP position, association significance, and SNP tier.
- Allelic heterogeneity: Secondary SNP enrichment in CRE data sets was tested against background and primary-tier SNPs using the SNP-tier indicator variable.Enrichment relative to background was assessed through the difference between β3 estimates.
- Insulator analysis: CTCF enrichment between independently associated SNPs for the same expression trait was evaluated against randomly selected cis-linked SNP pairs while controlling for genomic context.Controls included inter-SNP distance, position relative to the gene TSS, intervening recombination hotspots, and intervening TSSs.
Predicting replicating eQTLs with random forests
Random forests were used to predict whether eQTLs replicate within cell types, between cell types, or specifically to one cell type from genomic-location and CRE features. Ten-fold cross-validation measured generalization error, with ROC curves and AUC used for evaluation.
- Classifier design: Random forests predicted within-cell-type, between-cell-type, and cell-type-specific replication from eQTL genomic-location and cell-type CRE features.The ensemble approach was selected because it can capture interactions among features; performance was evaluated using 10-fold cross-validation and AUC.
Additional statistical analyses
The analyses used Fisher’s exact test, Wilcoxon’s rank sum test, and McNemar’s test for different categorical and interval-data comparisons.
- Statistical tests: Fisher’s exact, Wilcoxon’s rank sum, and McNemar’s tests were used for categorical, paired interval, and paired categorical data, respectively.The test choice depended on the structure of each comparison.