Interaction molecular QTL mapping discovers cellular and environmental modifiers of genetic regulatory effects

Abstract

Bulk-tissue molecular quantitative trait loci (QTLs) have been the starting point for interpreting disease-associated variants, and context-specific QTLs show particular relevance for disease. Here, we present the results of mapping interaction QTLs (iQTLs) for cell type, age, and other phenotypic variables in multi-omic, longitudinal data from the blood of individuals of diverse ancestries. By modeling the interaction between genotype and estimated cell-type proportions, we demonstrate that cell-type iQTLs could be considered as proxies for cell-type-specific QTL effects, particularly for the most abundant cell type in the tissue. The interpretation of age iQTLs, however, warrants caution because the moderation effect of age on the genotype and molecular phenotype association could be mediated by changes in cell-type composition. Finally, we show that cell-type iQTLs contribute to cell-type-specific enrichment of diseases that, in combination with additional functional data, could guide future functional studies. Overall, this study highlights the use of iQTLs to gain insights into the context specificity of regulatory effects.

Keywords: DNA methylation; cell-type composition; gene expression; gene-environment interaction; interaction QTL.

Figure 1

Study design and overview of the estimated cell-type proportions (A) Illustration of the study design and data types profiled for 1,319 individuals. (B) Graphical illustration of cell-type deconvolution. (C) Correlation of cell-type proportions. Exam 5 data from three sources were used: data estimated with the CIBERSORT method from PBMC gene expression, data estimated with the Houseman method from whole-blood DNA methylation, and cell counts measured by flow cytometry. (D) Sources of variability in estimated cell-type proportions with CIBERSORT and the Houseman method, gene expression from PBMCs, and DNA methylation from whole-blood according to exam 5 data. The median of the total explained variation is calculated across all the tested cell types, genes, and CpG sites. A gray dashed line denotes 1% of the total explained variance. Error bars denote the lower and upper quartile of the total explained variation.

Figure 2

Discovery of cell-type ieQTLs and imeQTLs (A) Illustration of the approach used for mapping cell-type interaction molQTLs in MESA. (B) Number of significant cell-type ieQTLs and imeQTLs combined across exams (FDR < 0.05 in exam 1 or exam 5 data) stratified by direction of the iQTL effect. (C) Reproducibility of cell-type iQTLs with a positive or negative direction of effect; one of the exams was used for discovery and the other for validation, and vice versa. The proportion of true positives (π₁ statistic) is used as a measure of reproducibility. (D) Sharing among cell-type ieQTLs and cell-type imeQTLs with a positive or negative direction of effect on the basis of exam 5 data is quantified as the proportion of true positives (π₁). The size of the square represents the correlation between the two estimated cell-type proportions measured via the absolute value of the Pearson correlation coefficient (r). (E) Sharing between CD4 T cell imeQTLs (query set) and cell-type ieQTLs (validation set) combined across exams is quantified as the proportion of CD4 T cell imeQTLs that have a positive direction of effect and are in LD (r² ≥ 0.5) with ieQTLs that have either positive or negative direction of effect from the given validation set. The p value shows the significance of the odds that CD4 T cell imeQTL with a positive direction of effect overlaps with a cell-type ieQTL with a positive direction of effect as compared to the odds that a CD4 T cell imeQTL with a positive direction of effect overlaps with a cell-type ieQTL with a negative direction of effect. (F) Example of a cell-type iQTL (rs774358) affecting both the expression levels of a gene (C9orf72) and a nearby CpG site (cg01126010). The p value of the interaction effect from the linear model fitted with TensorQTL is shown.

Figure 3

Replication and functional enrichment analysis of cell-type iQTLs (A) Replication of ieQTLs with a positive direction of effect in eQTL datasets from purified cell types from the eQTL Catalogue was based on effect size in eQTL data and allelic concordance. Highlighted are up to five datasets with absolute median effect size (beta) > 0.15 in the eQTL dataset and the proportion of QTLs with the same allelic direction >0.75 for B cell ieQTLs or >0.8 for other cell-type ieQTLs. Numerical results for all reference cell types are reported in Table S2. (B) Functional enrichment analysis performed with GoShifter shows the delta overlap, which is the difference between the observed proportion of loci overlapping a cCRE and the null for cell-type ieQTLs (upper panel) overlapping cCRE with high H3K27ac and for cell-type imeQTLs (lower panel) overlapping cCRE-dELSs. Negative delta overlap denotes that a smaller proportion of iQTL variants overlap with a cCRE than in the null distribution. Error bars denote the lower and upper quartile of the delta overlap. ^∗∗∗Significant association (adjusted p < 0.05) after correction for the number of target cell types with cCRE data, the number of cell types tested for interaction effect, and the number of groups of direction of effect. Numerical results for all reference cell types are reported in Table S3.

Figure 4

Trait iQTLs and mediated moderation (A) Number of significant trait ieQTLs and imeQTLs in exam 1 and exam 5 (FDR < 0.25) by direction of the iQTL effect. For numeric traits, direction of the iQTL effect is defined as (1) positive—genotype effect size increases depending on the trait or (2) negative—genotype effect size decreases depending on the trait. For binary traits, the direction of the effect is defined as (1) no effect in one—nominally non-significant genotype effect in one of the groups, (2) magnitude difference—nominally significant genotype effect in both groups with the same sign of the estimate, or (3) opposite effect—nominally significant genotype effect in both groups with the opposite sign of the estimate. (B) Example of smoking-current ieQTL for AHRR (upper plot) and age imeQTL for cg06953865 (lower plot); the p value of the interaction effect from the linear model fitted with TensorQTL is shown. (C) Inflation of G × monocyte effect among age ieQTLs and G × neutrophil effect among age imeQTLs in exam 5 data by direction of age iQTL effect. λ is the inflation factor. (D) Schema of the mediated-moderation approach, where the moderation effect of age on the genotype to DNAm association is mediated by changes in neutrophil proportions. The mediated-moderation effect is described by the G × age → G × neutrophil → DNAm path. p value histogram of the average causal mediation effect (ACME) of the G × neutrophil effect mediating the G × age effect on DNAm for 32 age imeQTLs with a positive or negative direction of effect.

Figure 5

Cell-type interaction QTLs and relevance for diseases (A) Relevance of cell-type ieQTLs (FDR < 0.25) and cell-type imeQTLs (FDR < 0.05) for selected cardiometabolic and immune diseases compared to height. For each of the cell-type iQTLs, we calculated the odds ratio (OR) as the ratio of the odds that an iQTL would colocalize with a cardiometabolic or immune disease to the odds that an iQTL would colocalize with height. For testing the significance of the OR, at least 10 loci had to be tested for colocalization; otherwise, the significance is noted as NA (not available). Bonferroni correction was applied separately for cell-type ieQTLs and cell-type imeQTLs. NS: not significant. (B) Colocalization between GWAS for RA by Okada et al., NK-cell ieQTLs for SYNGR1, and imeQTLs for a nearby CpG site cg19713460, shown as regional association plots. The highlighted region is depicted at the top and shows the location of rs909685, the lead GWAS variant for RA, and the CpG site relative to SYNGR1. (C) Association plot for the NK cell ieQTL for SYNGR1 and the NK-cell imeQTL for cg19713460. The p value of the interaction effect from the linear model fitted with TensorQTL is shown. Dots are colored on the basis of the genotype of rs909685. Data in (B) and (C) are from exam 1, where we observed the lowest interaction p values.