
==== Front
HGG Adv
HGG Adv
Human Genetics and Genomics Advances
2666-2477
Elsevier

S2666-2477(24)00088-5
10.1016/j.xhgg.2024.100348
100348
Article
Extensive co-regulation of neighboring genes complicates the use of eQTLs in target gene prioritization
Tambets Ralf 1
Kolde Anastassia 23
Kolberg Peep 1
Love Michael I. 45
Alasoo Kaur kaur.alasoo@ut.ee
16∗
1 Institute of Computer Science, University of Tartu, Tartu, Estonia
2 Institute of Genomics, University of Tartu, Tartu, Estonia
3 Institute of Mathematics and Statistics, University of Tartu, Tartu, Estonia
4 Department of Biostatistics, University of North Carolina at Chapel Hill, Chapel Hill, NC, USA
5 Department of Genetics, University of North Carolina at Chapel Hill, Chapel Hill, NC, USA
∗ Corresponding author kaur.alasoo@ut.ee
6 Lead contact

29 8 2024
10 10 2024
29 8 2024
5 4 1003482 2 2024
27 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by/4.0/ This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
Summary

Identifying causal genes underlying genome-wide association studies (GWASs) is a fundamental problem in human genetics. Although colocalization with gene expression quantitative trait loci (eQTLs) is often used to prioritize GWAS target genes, systematic benchmarking has been limited due to unavailability of large ground truth datasets. Here, we re-analyzed plasma protein QTL data from 3,301 individuals of the INTERVAL cohort together with 131 eQTL Catalog datasets. Focusing on variants located within or close to the affected protein identified 793 proteins with at least one cis-pQTL where we could assume that the most likely causal gene was the gene coding for the protein. We then benchmarked the ability of cis-eQTLs to recover these causal genes by comparing three Bayesian colocalization methods (coloc.susie, coloc.abf, and CLPP) and five Mendelian randomization (MR) approaches (three varieties of inverse-variance weighted MR, MR-RAPS, and MRLocus). We found that assigning fine-mapped pQTLs to their closest protein coding genes outperformed all colocalization methods regarding both precision (71.9%) and recall (76.9%). Furthermore, the colocalization method with the highest recall (coloc.susie - 46.3%) also had the lowest precision (45.1%). Combining evidence from multiple conditionally distinct colocalizing QTLs with MR increased precision to 81%, but this was accompanied by a large reduction in recall to 7.1%. Furthermore, the choice of the MR method greatly affected performance, with the standard inverse-variance-weighted MR often producing many false positives. Our results highlight that linking GWAS variants to target genes remains challenging with eQTL evidence alone, and prioritizing novel targets requires triangulation of evidence from multiple sources.

Linking non-coding genetic variants from association studies to their causal target genes is an important challenge in human genetics. We evaluate three colocalization and three Mendelian randomization methods on a gold standard dataset of cis protein quantitative trait loci. We find that co-regulation between neighboring genes complicates causal gene prioritization.

Keywords

GWAS
Mendelian randomisation
colocalisation
eQTL
pQTL
==== Body
pmcIntroduction

Linking non-coding regulatory variants from genome-wide association studies (GWASs) to their causal target genes is a fundamental problem in human genetics. Several strategies have been developed to address this problem. First, variants can simply be assigned to their closest protein coding genes. Second, colocalization with gene expression quantitative trait loci (eQTLs) can be used to ensure that a single GWAS variant also regulates gene expression.1,2 Finally, Mendelian randomization (MR) can be used to assess if multiple conditionally distinct variants have proportional effects on gene expression and the GWAS trait.3,4,5 However, systematic comparison of these strategies has been limited by methodological differences between studies and lack of comprehensive ground truth datasets linking trait-associated genetic variants to their causal genes.

Even for the simple strategy of assigning each GWAS variant to the closest gene, different studies have yielded varying estimates of precision and recall depending on which gene-variant pairs are being used as the truth set.6,7,8,9 Among others, these studies include the locus-2-gene (L2G) model,6 the activity-by-contact (ABC) model,7 ProGeM,8 and polygenic priority score (PoPS).9 The L2G model used 445 “gold-standard-positive” genes selected manually based on domain knowledge and literature review.6 The ABC model was evaluated on an enhancer perturbation dataset consisting of 109 regulatory connections inferred from experimental data.7 The PoPS model used fine-mapped missense variants to define their ground truth set.9 Finally, ProGeM used two different ground truth datasets: (1) 227 metabolite GWAS hits each assigned to high-confidence causal genes based on literature evidence, and (2) 562 cis protein quantitative trait loci (cis-pQTLs) data from the INTERVAL10 study, assuming that the most likely causal gene responsible for each cis-pQTL signal was the gene coding for the protein.8 As expected, the closest gene approach produced different results in the four studies: 37% recall and 47% precision in the ABC study; 55% recall and 56% precision in the L2G study; and 48% recall and 46% precision in the PoPS study. In contrast, ProGeM achieved 76% precision for metabolite GWAS hits and 69% precision for cis-pQTLs.

Here, we used the same INTERVAL cis-pQTL ground truth dataset employed by ProGeM but expanded the analysis in multiple ways. First, we fine-mapped the cis-pQTL signals, allowing us to consider multiple conditionally distinct causal variants for each protein in the same cis region. Second, we used colocalization instead of simple eQTL lookup to link cis-pQTLs to putative target genes.1 Third, fine-mapped eQTL data from the eQTL Catalogue allowed us to perform colocalization at the resolution of individual signals instead of genomic regions.2,11 Finally, identifying gene-protein pairs with two or more colocalizing signals allowed us to evaluate five MR approaches for causal gene prioritization. The closest gene approach (72% precision) outperformed all colocalization methods. Combining colocalization with MR and restricting analysis to gene-protein pairs with two or more shared signals increased the precision of target gene identification from 45% to 81%, but this came with a large decrease in recall (from 46% to 7%). The reduction in recall was primarily driven by the small sample size of eQTL datasets that limited the power to detect secondary eQTL signals. Importantly, we found that cis-eQTLs often violated one or more MR assumptions and using robust inference methods that accounted for these violations was essential to avoid false positives.

Results

Overview of the experimental design

To compare different colocalization methods on real-world data, we integrated cis-eQTL data from 131 fine-mapped tissue-specific datasets from the eQTL Catalogue (n = 65–702 individuals, Table S1) with fine-mapped plasma protein QTL data from the INTERVAL cohort (n = 3,301).10,12 Since the cell type and/or tissue source of any given plasma protein is often unclear,13 we decided to perform the colocalization analysis in a tissue-agnostic manner. We used three Bayesian colocalization methods: coloc.abf,1 which assumes a single causal variant per locus; coloc.susie,2 which supports multiple fine-mapped causal variants; and colocalization posterior probability (CLPP)11 defined at the variant level. We considered a colocalization signal to be significant in a locus if the posterior probability of colocalization (PP4) was greater than 0.8 for coloc.abf or coloc.susie. For CLPP, we used the commonly used threshold of 0.1 as the CLPP value is not directly comparable to PP4 from coloc.abf and coloc.susie.9,14

To illustrate how the three colocalization methods work, we looked at the colocalization between periostin (POSTN [MIM: 608777]) gene expression in the GTEx fibroblast dataset (QTD000216, n = 483) and POSTN plasma protein abundance in the INTERVAL dataset. The colocalization was clearly detected by coloc.abf (PP4 = 0.975) (Figure 1B). Interestingly, at this locus coloc.susie detected a colocalization for two independent fine-mapped signal pairs. The first eQTL signal colocalized with the second pQTL signal (PP4 = 0.979) (Figure 1C) and the second eQTL signal colocalized with the fifth pQTL signal (PP4 = 0.944) (Figure 1D). The two independent signals can also be seen on the coloc.abf plot (Figure 1B), where they had proportional effects on gene and protein levels and thus did not interfere with the single causal variant assumption of coloc.abf. CLPP did not detect a colocalization at this locus, because for the first fine-mapped signal pair, the CLPP value was below the 0.1 threshold (Figure 1C) and for the second signal pair, SuSiE did not detect a credible set for the fifth pQTL signal.Figure 1 Overview of datasets and analysis methods

(A) Fine-mapped cis-eQTLs from 131 distinct datasets and cis-pQTLs from the INTERVAL cohort were retrieved from eQTL Catalogue release 6. The coloc.abf, coloc.susie, and CLPP methods were used to identify colocalizing cis-eQTLs and cis-pQTLs pairs.

(B) coloc.abf colocalization between POSTN gene expression in GTEx fibroblasts (QTD000216, n = 483) and POSTN protein abundance in plasma.

(C and D) coloc.susie colocalization for two pairs of fine-mapped eQTL-pQTL signals in the same datasets. The red dot represents the shared lead QTL variant for the first eQTL signal and the second pQTL signal. The two purple dots represent the strongly linked (r2 = 0.995) lead variants for the second eQTL signal and the fifth pQTL signal. Instead of marginal association summary statistics, coloc.susie uses log Bayes factors (LBFs) for each fine-mapped signal.

(E) MR with the pair of colocalizing eQTLs from (C) and (D) is consistent with a causal effect from gene expression to protein abundance (graph created by MRLocus). The error bars represent standard errors.

(F) All five MR methods produce similar effect estimates and confidence intervals for the example in (E). ∗MRLocus provides a Bayesian credible interval instead of a confidence interval.

For the 96 gene-protein pairs for which coloc.susie detected multiple independent colocalizing signals in a single dataset, we further tested five MR methods (Figures 1E and 1F) to check if the effect sizes of the distinct colocalizing genetic signals were consistent with a putative causal effect of gene expression on protein abundance. We started with inverse-variance-weighted Mendelian randomization (IVW-MR) implemented in the MendelianRandomization R package15 and the IVW-MR with delta weights16 to account for the uncertainty of instrument effects on exposure. To account for potential violations of MR assumptions in the eQTL data, we also included three other methods that include additional modeling considerations: multiplicative random-effects IVW-MR that models overdispersion heterogeneity between instruments,16,17 MRLocus that models the dispersion of instrument’s effects via the allelic spread parameter,3 and MR-RAPS that models overdispersion heterogeneity while also accounting for outlier instruments.18 In the case of POSTN, all MR methods yielded very similar causal effect estimates and confidence intervals (Figure 1F).

Evaluating colocalization for causal gene identification

We first counted the number of proteins that were found to colocalize with any gene in any of the eQTL datasets. Since coloc.abf is only able to work with the strongest signal in a locus, we restricted coloc.susie to use only the first pQTL signal and CLPP to use only the first pQTL credible set for each protein for an objective comparison of the three methods. Even with these restrictions, coloc.susie found the largest number of proteins to colocalize (482 of the 793 tested), 57 of which were not discovered using the other two methods (Figure 2A). CLPP, on the other hand, found just 183 colocalizations, all of which were also detected by the other methods. Considering all independent pQTL signals increased the advantage of coloc.susie even further (Figure S1).Figure 2 Comparing the performance of three colocalization methods in causal gene identification

(A) The histogram on the left shows the number of proteins with at least one colocalizing eQTL detected by the CLPP, coloc.abf, and coloc.susie methods. The histogram on the right is the UpSetR19 plot showing the overlap of the colocalization events detected by the three methods. For CLPP and coloc.susie, only the first fine-mapped pQTL signal was included in the analysis.

(B) The precision and recall of the three colocalization methods in causal gene identification relative to the closest gene approach.

(C) Overview of the independent fine-mapped pQTL credible sets detected for each protein. The top histogram shows the number of proteins with 1, 2, 3, 4, and 5 or more credible sets detected. The density plots at the bottom show the distance from the fine-mapped pQTL credible set lead variant to the gene body of the corresponding protein coding gene.

(D) The precision and recall of the closest gene and coloc.susie methods as a function of the credible set index. CLPP, colocalization posterior probability.

Relying on the central dogma, we considered only significant colocalization signals between a cis-pQTL and the gene coding for the protein to be true positives (TPs) and all other significant signals to be false positives (FPs). Since our analysis was limited to cis-eQTLs, all detected false positive genes had to be located at most +/− 1 Mb from the cis-pQTL lead variant. However, on average, false positive genes were located farther away from the cis-pQTL lead variant than true positive genes (Figure S2). If a protein had at least one significant pQTL in the dataset but a colocalizing signal between it and its coding gene was not found, this was considered a false negative (FN). This allowed us to assess the recall (TP/(TP + FN)), the percentage of analyzed proteins found to colocalize with the gene coding for them and the precision (TP/(TP + FP)), the percentage of correct protein-gene pairs among all unique protein-gene pairs of each method.

Using only the first signal for each protein, coloc.susie found the correct gene for 367 of the 793 proteins with purity-filtered credible sets (46.3% recall). The recall was similar for signals located within the gene body (45.3%) or outside (47.9%). Coloc.abf was able to match a similar percentage of them to the coding gene (44.1%) at a slightly higher precision (48.1% vs. 45.1%) (Figure 2B). The variant-based approach of CLPP was the most precise of the three (68.5%), but yielded the correct gene for less than a fifth of the proteins (17.5%) (Figure 2B). However, all three colocalization methods were outperformed by a simple heuristic that assigned each pQTL to the gene body of the closest protein coding gene (76.9% recall, 71.9% precision) (Figure 2B). Notably, using distance to the closest transcription start site (TSS) instead of gene body decreased both precision (67.9%) and recall (68.1%), suggesting that some pQTLs might alter protein abundance in a TSS-independent manner (e.g., missense or 3′ UTR variants altering protein or mRNA stability). We did not find strong evidence that in the case of plasma pQTLs, restricting colocalization to specific cell types or tissues could be used to increase precision without significantly reducing recall (Note S1).

We then speculated that the good performance of the closest gene approach could be caused by the strongest pQTLs being located close to or within their target genes. Indeed, the median distance from the first fine-mapped pQTL signal to the corresponding protein coding gene was 0 base pairs (bp), meaning that most primary pQTLs were located within the gene body of the corresponding gene. This increased to 11,084 bp for the fifth and further signals with a wide spread (Figure 2C). We observed that both the precision (71.9% vs. 58.5%) and recall (76.9% vs. 62.0%) of the closest gene approach decreased slightly for tertiary and further pQTL signals (Figure 2D). The precision of the coloc.susie method was less affected by the pQTL signal index, but still always remained below the closest gene approach (e.g., 42.2% vs. 58.5% for third signals). Furthermore, recall of the colocalization approach decreased significantly for secondary pQTL signals (Figure 2D), suggesting that there might be less power to detect colocalizations at secondary signals due to their smaller effect sizes.

Protein QTLs detected on the SomaLogic platform are known to be susceptible to aptamer binding artifacts whereby missense variants in the protein sequence might alter aptamer binding affinity without changing protein abundance, thus giving rise to false positive pQTLs.10,13,20,21 We expected these variants to reduce the recall of our colocalization approach without substantially affecting precision as artifactual pQTLs should be less likely to colocalize with eQTLs. To test this, we restricted our analysis to 236 confidently fine mapped primary pQTL signals (PIP >0.8), 61 of which were missense variants and 175 were not. Only 21 of 61 missense variants colocalized with at least one eQTL (47.1% precision, 21.3% recall). In contrast, 95 of the 175 non-missense variants colocalized with at least one eQTL (55.3% precision, 44.6% recall). Interestingly, 13 of 21 colocalizing missense variants also overlapped a LeafCutter splicing QTL credible set from the eQTL Catalogue whereas only three of 40 non-colocalizing missense variants did. This is consistent with reports that some missense variants might also disrupt RNA splicing,22 thus potentially giving rise to weak but detectable eQTL signals.23 Of note, such sQTLs could still induce aptamer binding artifacts, making the sQTL-pQTL overlaps difficult to interpret.24

Evidence from multiple colocalizing eQTLs improves precision

Fine-mapping allowed us to identify multiple conditionally distinct cis-QTLs for both proteins and genes. In the INTERVAL dataset, we detect two or more pQTLs for 445 (56%) proteins with 71 (9%) having five or more independent signals in the cis region (Figure 3A). In the eQTL Catalog, the number of genes with multiple independent eQTLs depended on the sample size, but even in the group of datasets with the largest sample size (n > 350) only 20% of the genes had multiple independent fine-mapped eQTLs (Figure 3B). This suggests that most current eQTL datasets are too small to effectively fine-map multiple independent signals. This is consistent with previous analyses conducted by the GTEx, MetaBrain, and AdipoExpress projects, where the number of secondary eQTLs was strongly dependent on the eQTL sample sizes and reached more than 50% for tissues with the largest sample size.25,26,27Figure 3 Using Mendelian randomization to assess effect size concordance between colocalizing signal pairs

(A) Number of proteins with evidence of multiple signals (as indicated by credible set [CS] index) in the INTERVAL study (n = 3,301).

(B) Average number of genes with evidence of multiple signals (as indicated by CS index), stratified by the eQTL Catalogue dataset sample size.

(C) Two colocalizing QTL signals between glutathione S-transferase, omega-2 (GSTO2 [MIM: 612314]) gene expression in the liver (QTD000266, n = 208) and glutathione S-transferase, omega-1 (GSTO1 [MIM: 605482]) protein abundance in plasma. Standard IVW-MR detects a highly significant effect (p value <3∗10−290) that is not supported by the data.

(D) The signal from (C) analyzed with MRLocus. MRLocus 80% credible interval overlaps zero with wide allelic spread, because the two QTLs have inconsistent effects on gene expression and protein abundance.

(E) Significant MR signal between selectin L (SELL [MIM: 153240]) gene expression in blood (QTD000549, n = 195) and SELL protein abundance in plasma.

(F) Significant MR signal between selectin P (SELP [MIM: 173610]) gene expression in BLUEPRINT neutrophils (QTD000026, n = 196) and SELL protein abundance in plasma.

The error bars on panels C–F represent standard errors.

Coloc.susie is able to consider all pairwise colocalizations between independent fine-mapped QTLs in the same cis region (Figures 1C and 1D). Across all eQTL Catalogue datasets, we detected 321 gene-protein pairs with two or more colocalizing QTLs in the same cis region. The same gene-protein pair was often detected in multiple eQTL datasets, with 96 of them being unique (Table 1). These low numbers are primarily caused by a lack of secondary fine-mapped eQTL signals detected in the eQTL Catalogue (Figure 3B). As a result, this approach had a recall of only 8.6%, but precision increased from 45.1% (Figure 2B) to 70.8% (Table 1), suggesting that observing multiple colocalizing signals can significantly improve the precision of target gene identification.Table 1 Comparison of the MR methods used

Method	Precision	Recall	TP with positive slope	FP with positive slope	
coloc.susie (multiple signals), no MR	68/96 (70.8%)	68/793 (8.6%)	N/A	N/A	
IVW-MR	65/92 (70.7%)	65/793 (8.2%)	240/263 (91.3%)	25/48 (52.1%)	
IVW-MR (delta weights)	63/88 (71.6%)	63/793 (7.9%)	232/252 (92.1%)	23/44 (47.7%)	
IVW-MR (delta weights, random model)	56/69 (81.1%)	56/793 (7.1%)	205/215 (95.3%)	10/24 (41.7%)	
MRLocus	48/63 (76.2%)	48/793 (6.1%)	196/203 (96.6%)	10/21 (47.6%)	
MR-RAPS	64/87 (73.6%)	64/793 (8.1%)	229/254 (90.2%)	22/43 (51.2%)	
TP with positive slope - the fraction of true positive gene-protein-dataset triplets for which MR fitted a positive slope. FP with positive slope - the fraction of false-positive gene-protein-dataset triplets for which MR fitted a positive slope. While precision and recall are calculated at the level of gene-protein pairs, the direction of the MR slope is estimated separately in each gene-protein-dataset triplet. The methods with the highest precision and the highest proportion of true positives with a positive slope are shown in bold. We used the 80% credible interval for MRLocus and the 80% confidence interval for all other methods. N/A, not applicable.

Effect size concordance between multiple colocalizing QTLs

When focusing on gene-protein pairs with two or more independent colocalizing QTL signals, we noticed that the effect sizes of these QTLs were often discordant (Figures 3C and S3). We hypothesized that excluding these strongly discordant gene-protein pairs that are incompatible with a causal effect from gene expression to protein abundance might further increase precision. We first used inverse-variance weighted Mendelian randomization (IVW-MR) to identify these discordant pairs. Unexpectedly, we found that filtering colocalization results for significant MR slope (80% confidence interval not intersecting zero) had only minimal effect on identified gene-protein pairs (92 of 96 remained) and virtually no effect on precision and recall (Figure 3C; Table 1). These results remained robust to using more stringent filtering criteria (e.g., all 92 pairs remained at 95% confidence interval [Table S2]). Using IVW-MR with delta weights to better model large standard errors on the exposure typical for eQTL data also had minimal effect (Table 1). When consulting the literature, we realized that this behavior is likely driven by an assumption of the standard IVW-MR model that all instruments included in the analysis provide a (noisy) estimate of the same underlying causal effect. Violations of this assumption (Figures 3C and S3) can lead to underestimation of standard errors (Figure 3C) and overestimation of statistical significance.16,18,28

We tested three robust MR methods that explicitly model overdispersion heterogeneity or dispersion of instrument’s effects: multiplicative random-effects IVW-MR,16 MRLocus,3 and MR-RAPS.18 For the GSTO2-GSTO1 example (Figure 3C), both the random-effects IVW-MR and MRLocus now correctly inferred a null effect while MR-RAPS still detected a significant effect (Figures 3D and S4). The same was true for several other examples (Figure S5). Overall, random-effect IVW-MR (81.1% precision, 7.1% recall) performed slightly better than MRLocus (76.2% precision, 6.1% recall) while MR-RAPS (73.6% precision, 8.1% recall) performed similarly to the standard IVW-MR method (Table 1). On closer inspection, the poor performance of MR-RAPS seemed to be caused by its outlier detection feature that often excluded one of the two instruments from the analysis (Figures S4 and S5).

Concordance in gene and protein effect size direction

An assumption that we can make when working with gene-protein pairs is that for true causal relationships, the variants that increase gene expression should also increase protein abundance (i.e., have a positive MR slope). We split the gene-protein pairs identified by the five MR methods into TPs if the gene coded for the protein and FPs otherwise. We found that MRLocus had the highest fraction of true positive pairs with a positive slope (96.6%) closely followed by random-effect IVW-MR (95.3%) (Table 1). The other three methods had lower effect size direction concordance that ranged from 90.2% (MR-RAPS) to 92.1% (delta-weighted IVW-MR) (Table 1), indicating that explicit modeling of overdispersion heterogeneity not only increases precision, but the gene-protein pairs detected with robust methods (MRLocus and random-effect IVW-MR) also more often display the expected direction of effect. In contrast, ∼50% of the false positive gene-protein pairs from the five MR methods had a positive slope. In some cases, the remaining FPs from the MRLocus and random-effect IVW-MR analysis with positive slopes reflected genes in the same cis locus where two or more instruments had highly concordant effects on both genes, highlighting how strong local co-regulation can confuse even the best causal inference approaches (Figures 3E and 3F). Finally, although we had too few gene-protein pairs to stratify the analysis by cell types and tissues (Table S9), we did detect one example where contrasting MR slopes and allelic heterogeneity estimates between cell types helped to prioritize the likely causal cell type (Figure S6).

Increasing the power of cis-MR via eQTL meta-analysis

Although MR with two or more colocalizing signals increased the precision of target gene identification to 81%, this came at the cost of significantly reduced recall (7.1%). The primary reason for detecting few multi-signal gene-protein pairs in our analysis is the relatively small sample size of eQTL datasets (n = 65–702) in the eQTL Catalogue. To test if the recall of the MR methods could be increased by meta-analyzing eQTLs across multiple studies from the same tissue, we obtained the cis-eQTL summary statistics from the AdipoExpress project, a meta-analysis of five subcutaneous adipose tissue studies (n = 2,344). Although AdipoExpress was not able to use SuSiE for fine-mapping due to the risk of FPs,29 they used all-but-one conditional analysis30 to identify conditionally distinct signals. After converting these conditional summary statistics to approximate Bayes factors (see subjects, material, and methods), we performed colocalization with INTERVAL pQTLs using the same workflow that we previously used for the eQTL Catalogue datasets.

We compared the colocalization that we detected in the AdipoExpress dataset with those from the best-powered adipose tissue dataset from the eQTL Catalogue (TwinsUK, n = 381). We found that the number of colocalizing gene-protein pairs increased by approximately 2-fold for both coloc.abf (from 98 to 195) and coloc.susie (104–242, first signal only). Consistent with the observation that larger sample sizes increase the power to detect secondary eQTL signals (Figure 3B),25,26,27 the number of multi-signal gene-protein pairs increased from eight to 35 (4.4-fold), but only one multi-signal pair was shared between the two analyses. Furthermore, 18 of 35 gene-protein pairs also had a significant MR effect in the AdipoExpress dataset (95% precision) (Table S3) and seven of eight had a significant MR effect in the TwinsUK dataset (87% precision) (Table S4). However, overall recall remained low (2.3% for AdipoExpress and 0.9% for TwinsUK, Tables S3 and S4), potentially because adipose tissue is unlikely to be the causal tissue for many plasma cis-pQTLs.

Discussion

A fundamental problem in human complex traits genetics is linking primarily non-coding GWAS hits to their causal target genes. Here, we used fine-mapped cis-pQTLs to systematically evaluate the performance of eQTL colocalization methods in identifying causal target genes. Our key assumption was that the causal gene responsible for a cis-pQTL signal should be the gene coding for the protein. Our results indicate that eQTL colocalization approaches, when performed systematically against very large eQTL databases such as the eQTL Catalogue, have generally low precision (∼50%) in identifying the correct target genes. This seems to be primarily driven by horizontal pleiotropy whereby the same eQTL variants are associated with the expression level of multiple genes located in the same cis locus (+/− 1 Mb). We also found that precision can be improved (up to 81%) when combining multiple colocalizing QTL signals in an MR framework to explicitly consider the concordance of the causal effect estimates provided by independent genetic variants. This agrees with other recent studies, affirming that combining colocalization with MR reduces confounding by linkage disequilibrium (LD) and improves the sensitivity and specificity of identifying biologically relevant targets.31,32,33

In our analysis, the closest gene approach had very high precision (71.9%), which is higher than typically seen in other studies that benchmark methods for causal gene prioritization in the GWAS setting (range 46%–56%).6,7,9 Part of the reason could be that primary pQTLs might be much closer to their target genes than typical GWAS hits are (Figure 1C), matching a similar observation for eQTL variants and GWAS hits.34 Indeed, we observed that the precision of the closest gene approach dropped to 58.5% for tertiary and further pQTL signals that were more often located outside the gene body (Figure 1D). This suggests that in realistic GWAS target prioritization applications, colocalization and closest gene approach might achieve similarly moderate precision of ∼50%. Second, our closest gene approach was based on the distance to the gene body as opposed to the closest TSS chosen by some other studies.6 Indeed, using TSS instead of gene body to define the closest genes decreased the precision to 67.8% (Figure 2B). While restricting colocalization to a small number of trait-relevant tissues or cell types is sometimes used to reduce FPs,35 we found that this can significantly reduce recall (Note S1). Thus, we would recommend using all available data for initial colocalization analysis followed by complementary methods to prioritize likely causal genes.

Our results highlight the challenges of using gene expression levels as exposures in MR. The primary concern is that eQTL variants often do not satisfy the exclusion restriction assumption of the MR framework,36 which states that genetic variants affect the outcome only through their effect on the gene expression level (exposure) included in the model. This assumption can be violated in at least two ways. First, we might be looking at the right gene but in the wrong context. The true causal effect of the gene expression on the outcome might be mediated in some other tissue, cell type, or context that was not included in the analysis. In this scenario, failure to model allelic spread or overdispersion heterogeneity may inflate the significance of the standard IVW-MR estimates.3,17 Second, due to horizontal pleiotropy, we might be looking at the wrong gene and the actual causal effect might be mediated by another gene for which the eQTL variants have highly correlated effects (e.g., SELL and SELP genes on Figures 3E and 3F). Thus, we caution against interpreting significant cis-MR slopes as direct evidence of the causal effect of gene expression levels on the outcome. Rather, we prefer to use MR to exclude exposures that are clearly inconsistent with a causal effect in the tested cell type or tissue. Finally, our results reinforce the need to include both positive and negative controls in MR analysis and use visualization approaches to assess model fit.37

A promising approach that we did not evaluate here is multivariable MR, which jointly models the expression levels of all nearby genes.5,17 However, multivariable MR requires that the number of genetic instruments included in the analysis equals or exceeds the number of exposures, which is unrealistic for cis-eQTLs from large compendia containing hundreds of cell types and tissues.38 Furthermore, multivariable MR can identify the correct causal gene only if the right gene in the right context (cell type or tissues), or a sufficient proxy context, is included in the model as one of the exposures. There also needs to be sufficient phenotypic heterogeneity between the different exposures included in the model (i.e., genetic variant effects vary between the different genes).17 Thus, multivariable MR is unlikely to completely resolve the exclusion restriction assumption violations that we have observed here.

Choosing the closest gene almost always outperformed eQTL colocalization when identifying causal genes responsible for cis-pQTLs. Furthermore, even though MR with multiple independent eQTLs did outperform the closest gene approach in terms of precision (81.2% vs. 71.9%), this came at the cost of a significant reduction in recall (76.9% vs. 7.1%). This reduction in recall was primarily driven by limited power to detect secondary eQTL signals in existing datasets (Figure 3B). Consistent with this hypothesis, we found that using adipose tissue cis-eQTL conditional meta-analysis summary statistics from the AdipoExpress project (n = 2,344) instead of TwinsUK (n = 415) increased the recall of eQTL cis-MR by 2.5-fold from 0.9% to 2.3% (Tables S3 and S4). Thus, a promising avenue to improve the recall of eQTL MR is to increase the sample sizes of eQTL datasets by either collecting new samples or performing meta-analysis across multiple existing datasets. A potential added benefit is that secondary eQTLs might represent more distal context-specific effects that are more likely to overlap disease GWAS hits.34 Finally, successful target gene prioritization will likely require triangulation of evidence from multiple genetic and non-genetic sources. Fortunately, multiple competing statistical models are currently actively being developed to support this integration (e.g., L2G6 and PoPS9).

Subjects, material, and methods

Datasets used in the analysis

We downloaded eQTL summary statistics and fine-mapping results for 34 studies from the eQTL Catalogue (release 6) FTP server (https://www.ebi.ac.uk/eqtl/Data_access/).26,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67 The genotype and protein abundance data from the INTERVAL cohort10 were downloaded from EGA (accessions EGA:EGAD00010001544 and EGA:EGAD00001004080) after access was approved by the “Plasma pQTLs in INTERVAL cohort” data access committee. The INTERVAL study comprises about 50,000 participants nested within a randomized trial of varying blood donation intervals.68 Between mid-2012 and mid-2014, blood donors aged 18 years and older were recruited at 25 centers of England’s National Health Service Blood and Transplant (NHSBT). All participants gave informed consent before joining the study and the National Research Ethics Service approved this study (11/EE/0538). Participants completed an online questionnaire including questions about demographic characteristics (for example, age, sex, ethnicity), anthropometry (height, weight), lifestyle (for example, alcohol and tobacco consumption), and diet. For SomaLogic assays, two non-overlapping subcohorts of 2,731 and 831 participants were randomly selected from INTERVAL. After genetic quality control, 3,301 participants remained for analysis.10

INTERVAL pQTL data processing

Genotype imputation

The imputed genotypes from the INTERVAL cohort were based on the GRCh37 coordinates, but all fine-mapping results from the eQTL Catalogue used GRCh38 coordinates. We initially used CrossMap.py69 to convert imputed genotypes to GRCh38 coordinates but found that this approach caused some artifactual fine-mapping results due to variants lost in the lift-over process. To avoid these issues, we extracted genotyped variant positions of Affymetrix Axiom UK Biobank array from the INTERVAL imputed genotype files and re-imputed genotypes to the 1000 Genomes 30x on GRCh38 reference panel with the eQTL-Catalogue/genimpute v23.07.1 workflow. The same workflow was previously used to impute genotypes for all eQTL Catalogue datasets that used genotyping microarrays.12

Protein data processing and association testing

We downloaded the pre-processed SomaLogic protein abundance data and aptamer metadata from EGA (EGAD00001004080). We applied inverse normal transformation to the protein abundance data and used the g:Profiler70 web tool to map SomaLogic protein names to Ensembl gene ids. The metadata for SomaLogic aptamers and their mapping to Ensembl gene ids can be downloaded from Zenodo (https://doi.org/10.5281/zenodo.7808390). The cis-pQTL analysis and fine mapping were conducted using the eQTL-Catalogue/qtlmap v23.02.1 workflow as described previously.12 In all downstream analyses, we used 793 proteins that had at least one purity-filtered SuSiE credible set.

AdipoExpress data processing

The AdipoExpress project performed cis-eQTL meta-analysis across five subcutaneous adipose tissue studies (total n = 2,344).27 We downloaded AdipoExpress summary statistics from https://mohlke.web.unc.edu/data/adipoexpress/. Since AdipoExpess analysis used the GRCh37 reference genome, we first converted the variant positions to GRCh38 coordinates with the MungeSumstats R package.71 Since using SuSiE to identify conditionally distinct signals is prone to false positives in a meta-analysis setting,29 AdipoExpress used all-but-one conditional meta-analysis and also released summary statistics for distinct signals conditioned on all other significant signals in the same cis region. To use these results with coloc.susie, we first converted the all-but-one conditional betas and standard errors to log approximate Bayes factors (LABFs) using the process.dataset() function from the coloc R package. Subsequently, we used the LABFs in place of the log Bayes factors in the coloc.susie method.

Colocalization between cis-eQTLs and cis-pQTLs

We downloaded fine-mapped cis-eQTL summary statistics for 131 datasets of the eQTL Catalogue release 6 from eQTL Catalogue FTP server.12 We ran colocalization analyses pairwise between all eQTL datasets and pQTL data from the INTERVAL study. We set the cis-window for each locus at 2 million base pairs centered at the TSS. The Nextflow workflow implementing the CLPP, coloc.susie, and coloc.abf colocalization methods is available from GitHub (https://github.com/ralf-tambets/coloc). The workflow assigned colocalization probabilities for each cis-pQTL locus in the INTERVAL dataset and each cis-eQTL locus of each eQTL dataset. We only included protein coding genes in the analysis as these are much more likely to be the causal genes for pQTLs. We also excluded all protein complexes from the INTERVAL dataset as their abundance could be influenced by all their constituents independently.72,73

CLPP

We calculated CLPP as described previously.14 Briefly, we joined the data from the eQTL study and the pQTL study by the variant name. We calculated CLPP for each variant by multiplying the posterior inclusion probabilities from both studies and summed the resulting values up for each credible set in the eQTL dataset. A signal was significant if CLPP exceeded 0.1. This approach yielded 2,278 unique colocalizing gene-protein-dataset triplets (Table S5).

coloc.abf

We analyzed each dataset chromosome by chromosome by running coloc.abf1 on summary statistics (beta, standard error, MAF) between each eQTL gene and a subset of the cis-pQTLs that fell in the cis-window, unless more than 90% of the variants in the pQTL gene fell outside the cis-window. The prior probabilities that an SNP is associated with either trait were set at 1 × 10−4 and the prior probability that an SNP is associated with both traits was set at 5 × 10−6. A signal was significant if PP4 exceeded 0.8. This approach yielded 4,890 unique colocalizing gene-protein-dataset triplets (Table S6).

coloc.susie

Data preparation for coloc.susie2 was similar to that of coloc.abf with the exception that the input data consisted of SuSiE log Bayes factors (LBFs) for all fine-mapped signals instead of marginal betas and standard errors. We used the coloc.bf_bf function to calculate the colocalization posterior probabilities, which we ran with the same prior probabilities of association as for coloc.abf. A signal was significant if PP4 exceeded 0.8. This approach yielded 8,501 unique colocalizing gene-protein-dataset triplets (Table S7).

Summary analysis

To determine the closest gene to a given protein for benchmarking purposes, we found the lead variant for each credible set based on Z score and calculated the distance from it to the start and end coordinates of each protein coding gene. If the lead variant fell within the gene body, the distance to the gene was set to zero. In cases of equal closest distances, all tied genes were considered closest, with at most one of them being a true positive.

Mendelian randomization

We tested five different MR methods: default mr_ivw() method from the MendelianRandomization R package version 0.9.015; the same function with the weights argument set to “delta”16; the same function with the weights argument set to “delta” and the model argument set to “random”; the mr_raps() function from the MR-RAPS R package version 0.4.118 with the over.dispersion argument set to TRUE and the loss.function argument set to “tukey”; and the fitSlope() function from the MRLocus R package version 0.0.263. We considered the signal from a gene-protein-dataset triplet as significant if the 80% confidence interval did not include 0 (using MendelianRandomization and MR-RAPS) or if the 80% credible interval did not include 0 (using MRLocus).

Data and code availability

All eQTL summary statistics, fine-mapping results, and log Bayes factors are available from the eQTL Catalogue FTP server (https://www.ebi.ac.uk/eqtl/). The pQTL summary statistics and fine-mapping results from the INTERVAL cohort have also been deposited to the eQTL Catalogue under the accession QTD000584. The accession numbers for the individual-level genotypes and protein abundances from the INTERVAL cohort are EGA: EGAD00001004080 and EGA: EGAD00010001544. The code used for colocalization and Mendelian randomization analyses is available at https://github.com/ralf-tambets/coloc.

Acknowledgments

We thank S. Kasela for her helpful comments on the manuscript. The colocalization and Mendelian randomization analyses were performed at the High-Performance Computing Center, University of Tartu. We thank INTERVAL study participants; staff at recruiting NHSBT blood donation centers; and the INTERVAL Study Coordination team, Operations Team (led by R. Houghton and C. Moore) and Data Management Team (led by M. Walker). K.A., R.T., and P.K. were supported by the 10.13039/501100002301 Estonian Research Council (grant no. PSG415 ).

Author contributions

R.T. performed all colocalization and Mendelian randomization analyses presented in the paper. A.K. prepared the INTERVAL proteomics dataset for pQTL analysis. P.K. developed the genotype imputation workflow for low-coverage whole genome sequencing data. R.T., K.A., and M.I.L. interpreted the Mendelian randomization results. K.A. and R.T wrote the manuscript with input from all authors.

Declaration of interests

The authors declare no competing interests.

Web resources

Coloc, https://chr1swallace.github.io/coloc/index.html.

MRLocus, https://thelovelab.github.io/mrlocus/.

MendelianRandomization, https://github.com/cran/MendelianRandomization.

MR.RAPS, https://github.com/qingyuanzhao/mr.raps.

AdipoExpress, https://mohlke.web.unc.edu/data/adipoexpress/.

Colocalization workflow, https://github.com/ralf-tambets/coloc.

eQTL Catalogue FTP server, https://www.ebi.ac.uk/eqtl/.

eQTL-Catalogue/genimpute workflow, https://github.com/eQTL-Catalogue/genimpute.

eQTL-Catalogue/qtlmap workflow, https://github.com/eQTL-Catalogue/qtlmap.

Online Mendelian Inheritance in Man, https://omim.org/.

Supplemental information

Document S1. Figures S1–S6, Tables S2–S4, and Note S1

Data S1. Tables S1 and S5–S9

Document S2. Article plus supplemental information

Supplemental information can be found online at https://doi.org/10.1016/j.xhgg.2024.100348.
==== Refs
References

1 Giambartolomei C. Vukcevic D. Schadt E.E. Franke L. Hingorani A.D. Wallace C. Plagnol V. Bayesian test for colocalisation between pairs of genetic association studies using summary statistics PLoS Genet. 10 2014 e1004383
2 Wallace C. A more accurate method for colocalisation analysis allowing for multiple causal variants PLoS Genet. 17 2021 e1009440
3 Zhu A. Matoba N. Wilson E.P. Tapia A.L. Li Y. Ibrahim J.G. Stein J.L. Love M.I. MRLocus: Identifying causal genes mediating a trait through Bayesian estimation of allelic heterogeneity PLoS Genet. 17 2021 e1009455
4 van der Graaf A. Claringbould A. Rimbert A. BIOS ConsortiumWestra H.J. Li Y. Wijmenga C. Sanna S. Mendelian randomization while jointly modeling cis genetics identifies causal relationships between gene expression and lipids Nat. Commun. 11 2020 4930 5012 33004804
5 Porcu E. Rüeger S. Lepik K. eQTLGen ConsortiumBIOS ConsortiumSantoni F.A. Reymond A. Kutalik Z. Mendelian randomization integrating GWAS and eQTL data reveals genetic determinants of complex and clinical traits Nat. Commun. 10 2019 3300 31341166
6 Mountjoy E. Schmidt E.M. Carmona M. Schwartzentruber J. Peat G. Miranda A. Fumis L. Hayhurst J. Buniello A. Karim M.A. An open approach to systematically prioritize causal variants and genes at all published human GWAS trait-associated loci Nat. Genet. 53 2021 1527 1533 34711957
7 Fulco C.P. Nasser J. Jones T.R. Munson G. Bergman D.T. Subramanian V. Grossman S.R. Anyoha R. Doughty B.R. Patwardhan T.A. Activity-by-contact model of enhancer-promoter regulation from thousands of CRISPR perturbations Nat. Genet. 51 2019 1664 1669 31784727
8 Stacey D. Fauman E.B. Ziemek D. Sun B.B. Harshfield E.L. Wood A.M. Butterworth A.S. Suhre K. Paul D.S. ProGeM: a framework for the prioritization of candidate causal genes at molecular quantitative trait loci Nucleic Acids Res. 47 2019 e3 30239796
9 Weeks E.M. Ulirsch J.C. Cheng N.Y. Trippe B.L. Fine R.S. Miao J. Patwardhan T.A. Kanai M. Nasser J. Fulco C.P. Leveraging polygenic enrichments of gene features to predict genes underlying complex traits and diseases Nat. Genet. 55 2023 1267 1276 37443254
10 Sun B.B. Maranville J.C. Peters J.E. Stacey D. Staley J.R. Blackshaw J. Burgess S. Jiang T. Paige E. Surendran P. Genomic atlas of the human plasma proteome Nature 558 2018 73 79 29875488
11 Hormozdiari F. van de Bunt M. Segrè A.V. Li X. Joo J.W.J. Bilow M. Sul J.H. Sankararaman S. Pasaniuc B. Eskin E. Colocalization of GWAS and eQTL Signals Detects Target Genes Am. J. Hum. Genet. 99 2016 1245 1260 27866706
12 Kerimov N. Tambets R. Hayhurst J.D. Rahu I. Kolberg P. Raudvere U. Kuzmin I. Chowdhary A. Vija A. Teras H.J. eQTL Catalogue 2023: New datasets, X chromosome QTLs, and improved detection and visualisation of transcript-level QTLs PLoS Genet. 19 2023 e1010932
13 Pietzner M. Wheeler E. Carrasco-Zanini J. Cortes A. Koprulu M. Wörheide M.A. Oerton E. Cook J. Stewart I.D. Kerrison N.D. Mapping the proteo-genomic convergence of human diseases Science 374 2021 eabj1541
14 Kanai M. Ulirsch J.C. Karjalainen J. Kurki M. Karczewski K.J. Fauman E. Wang Q.S. Jacobs H. Aguet F. Ardlie K.G. Insights from complex trait fine-mapping across diverse populations Preprint at bioRxiv 2021 10.1101/2021.09.03.21262975
15 Yavorska O.O. Burgess S. MendelianRandomization: an R package for performing Mendelian randomization analyses using summarized data Int. J. Epidemiol. 46 2017 1734 1739 28398548
16 Burgess S. Bowden J. Integrating summarized data from multiple genetic variants in Mendelian randomization: bias and coverage properties of inverse-variance weighted methods Preprint at arXiv 2015 10.48550/arXiv.1512.04486
17 Patel A. Gill D. Shungin D. Mantzoros C.S. Knudsen L.B. Bowden J. Burgess S. Robust use of phenotypic heterogeneity at drug target genes for mechanistic insights: application of cis-multivariable Mendelian randomization toGLP1Rgene region Preprint at bioRxiv 2023 10.1101/2023.07.20.23292958
18 Zhao Q. Wang J. Hemani G. Bowden J. Small D.S. Statistical inference in two-sample summary-data Mendelian randomization using robust adjusted profile score Ann. Stat. 48 2020 1742 1769
19 Conway J.R. Lex A. Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties Bioinformatics 33 2017 2938 2940 28645171
20 Pietzner M. Wheeler E. Carrasco-Zanini J. Kerrison N.D. Oerton E. Koprulu M. Luan J. Hingorani A.D. Williams S.A. Wareham N.J. Langenberg C. Synergistic insights into human health from aptamer- and antibody-based proteomic profiling Nat. Commun. 12 2021 6822 6913 34819519
21 Ferkingstad E. Sulem P. Atlason B.A. Sveinbjornsson G. Magnusson M.I. Styrmisdottir E.L. Gunnarsdottir K. Helgason A. Oddsson A. Halldorsson B.V. Large-scale integration of the plasma proteome with genetics and disease Nat. Genet. 53 2021 1712 1721 34857953
22 Soemedi R. Cygan K.J. Rhine C.L. Wang J. Bulacan C. Yang J. Bayrak-Toydemir P. McDonald J. Fairbrother W.G. Pathogenic variants that alter protein code often disrupt splicing Nat. Genet. 49 2017 848 855 28416821
23 Kerimov N. Hayhurst J.D. Peikova K. Manning J.R. Walter P. Kolberg L. Samoviča M. Sakthivel M.P. Kuzmin I. Trevanion S.J. A compendium of uniformly processed human gene expression and splicing quantitative trait loci Nat. Genet. 53 2021 1290 1299 34493866
24 Tokolyi A. Persyn E. Nath A.P. Burnham K.L. Marten J. Vanderstichele T. Tardaguila M. Stacey D. Farr B. Iyer V. Genetic determinants of blood gene expression and splicing and their contribution to molecular phenotypes and health outcomes Preprint at medRxiv 2023 10.1101/2023.11.25.23299014
25 de Klein N. Tsai E.A. Vochteloo M. Baird D. Huang Y. Chen C.Y. van Dam S. Oelen R. Deelen P. Bakker O.B. Brain expression quantitative trait locus and network analyses reveal downstream effects and putative drivers for brain-related diseases Nat. Genet. 55 2023 377 388 36823318
26 GTEx Consortium The GTEx Consortium atlas of genetic regulatory effects across human tissues Science 369 2020 1318 1330 32913098
27 Brotman S.M. El-Sayed Moustafa J.S. Guan L. Broadaway K.A. Wang D. Jackson A.U. Welch R. Currin K.W. Tomlinson M. Vadlamudi S. Adipose tissue eQTL meta-analysis reveals the contribution of allelic heterogeneity to gene expression regulation and cardiometabolic traits Preprint at bioRxiv 2023 10.1101/2023.10.26.563798
28 Burgess S. Davey Smith G. Davies N.M. Dudbridge F. Gill D. Glymour M.M. Hartwig F.P. Kutalik Z. Holmes M.V. Minelli C. Guidelines for performing Mendelian randomization investigations: update for summer 2023 Wellcome Open Res. 4 2019 186 32760811
29 Kanai M. Elzur R. Zhou W. Global Biobank Meta-analysis InitiativeDaly M.J. Finucane H.K. Meta-analysis fine-mapping is often miscalibrated at single-variant resolution Cell Genom. 2 2022 100210
30 Brown M. Greenwood E. Zeng B. Powell J.E. Gibson G. Effect of All-but-One Conditional Analysis for eQTL Isolation in Peripheral Blood Genetics 223 2023 iyac162 10.1093/genetics/iyac162
31 Karim M.A. Ariano B. Schwartzentruber J. Roldan-Romero J.M. Mountjoy E. Hayhurst J. Buniello A. Mohammed E.S.E. Carmona M. Holmes M.V. Systematic disease-agnostic identification of therapeutically actionable targets using the genetics of human plasma proteins Preprint at medRxiv 2023 10.1101/2023.06.01.23290252
32 Hukku A. Sampson M.G. Luca F. Pique-Regi R. Wen X. Analyzing and reconciling colocalization and transcriptome-wide association studies from the perspective of inferential reproducibility Am. J. Hum. Genet. 109 2022 825 837 35523146
33 Zuber V. Grinberg N.F. Gill D. Manipur I. Slob E.A.W. Patel A. Wallace C. Burgess S. Combining evidence from Mendelian randomization and colocalization: Review and comparison of approaches Am. J. Hum. Genet. 109 2022 767 782 35452592
34 Mostafavi H. Spence J.P. Naqvi S. Pritchard J.K. Systematic differences in discovery of genetic effects on gene expression and complex traits Nat. Genet. 55 2023 1866 1875 37857933
35 Sobczyk M.K. Richardson T.G. Zuber V. Min J.L. Gaunt T.R. Paternoster L. eQTLGen Consortium, BIOS Consortium, GoDMC Triangulating molecular evidence to prioritize candidate causal genes at established atopic dermatitis loci J. Invest. Dermatol. 141 2021 2620 2629 10.1016/j.jid.2021.03.027 33901562
36 Davies N.M. Holmes M.V. Davey Smith G. Reading Mendelian randomisation studies: a guide, glossary, and checklist for clinicians BMJ 362 2018 k601 30002074
37 Hamilton F.W. Hughes D.A. Spiller W. Tilling K. Smith G.D. Non-linear mendelian randomization: evaluation of biases using negative controls with a focus on BMI and Vitamin D Preprint at bioRxiv 2023 10.1101/2023.08.21.23293658
38 Burgess S. Mason A.M. Grant A.J. Slob E.A.W. Gkatzionis A. Zuber V. Patel A. Tian H. Liu C. Haynes W.G. Using genetic association data to guide drug discovery and development: Review of methods and applications Am. J. Hum. Genet. 110 2023 195 214 36736292
39 Alasoo K. Rodrigues J. Mukhopadhyay S. Knights A.J. Mann A.L. Kundu K. HIPSCI ConsortiumHale C. Dougan G. Gaffney D.J. Shared genetic effects on chromatin and gene expression indicate a role for enhancer priming in immune response Nat. Genet. 50 2018 424 431 29379200
40 Chen L. Ge B. Casale F.P. Vasquez L. Kwan T. Garrido-Martín D. Watt S. Yan Y. Kundu K. Ecker S. Genetic Drivers of Epigenetic and Transcriptional Variation in Human Immune Cells Cell 167 2016 1398 1414.e24 27863251
41 Gutierrez-Arcelus M. Lappalainen T. Montgomery S.B. Buil A. Ongen H. Yurovsky A. Bryois J. Giger T. Romano L. Planchon A. Passive and active DNA methylation and the interplay with genetic variation in gene regulation Elife 2 2013 e00523
42 Lappalainen T. Sammeth M. Friedländer M.R. 't Hoen P.A.C. Monlong J. Rivas M.A. Gonzàlez-Porta M. Kurbatova N. Griebel T. Ferreira P.G. Transcriptome and genome sequencing uncovers functional variation in humans Nature 501 2013 506 511 24037378
43 Kilpinen H. Goncalves A. Leha A. Afzal V. Alasoo K. Ashford S. Bala S. Bensaddek D. Casale F.P. Culley O.J. Common genetic variation drives molecular heterogeneity in human iPSCs Nature 546 2017 370 375 28489815
44 Nédélec Y. Sanz J. Baharian G. Szpiech Z.A. Pacis A. Dumaine A. Grenier J.C. Freiman A. Sams A.J. Hebert S. Genetic Ancestry and Natural Selection Drive Population Differences in Immune Responses to Pathogens Cell 167 2016 657 669.e21 27768889
45 Quach H. Rotival M. Pothlichet J. Loh Y.H.E. Dannemann M. Zidane N. Laval G. Patin E. Harmant C. Lopez M. Genetic Adaptation and Neandertal Admixture Shaped the Immune System of Human Populations Cell 167 2016 643 656.e17 27768888
46 Schwartzentruber J. Foskolou S. Kilpinen H. Rodrigues J. Alasoo K. Knights A.J. Patel M. Goncalves A. Ferreira R. Benn C.L. Molecular and functional variation in iPSC-derived sensory neurons Nat. Genet. 50 2018 54 61 29229984
47 Buil A. Brown A.A. Lappalainen T. Viñuela A. Davies M.N. Zheng H.F. Richards J.B. Glass D. Small K.S. Durbin R. Gene-gene and gene-environment interactions detected by transcriptome sequence analysis in twins Nat. Genet. 47 2015 88 91 25436857
48 van de Bunt M. Manning Fox J.E. Dai X. Barrett A. Grey C. Li L. Bennett A.J. Johnson P.R. Rajotte R.V. Gaulton K.J. Transcript Expression Data from Human Islets Links Regulatory Signals from Genome-Wide Association Studies for Type 2 Diabetes and Glycemic Traits to Their Downstream Effectors PLoS Genet. 11 2015 e1005694
49 Schmiedel B.J. Singh D. Madrigal A. Valdovino-Gonzalez A.G. White B.M. Zapardiel-Gonzalo J. Ha B. Altay G. Greenbaum J.A. McVicker G. Impact of Genetic Polymorphisms on Human Immune Cell Gene Expression Cell 175 2018 1701 1715.e16 30449622
50 Jaffe A.E. Straub R.E. Shin J.H. Tao R. Gao Y. Collado-Torres L. Kam-Thong T. Xi H.S. Quan J. Chen Q. Developmental and genetic regulation of the human cortex transcriptome illuminate schizophrenia pathogenesis Nat. Neurosci. 21 2018 1117 1125 30050107
51 Ng B. White C.C. Klein H.U. Sieberts S.K. McCabe C. Patrick E. Xu J. Yu L. Gaiteri C. Bennett D.A. An xQTL map integrates the genetic architecture of the human brain’s transcriptome and epigenome Nat. Neurosci. 20 2017 1418 1426 28869584
52 Lepik K. Annilo T. Kukuškina V. eQTLGen ConsortiumKisand K. Kutalik Z. Peterson P. Peterson H. C-reactive protein upregulates the whole blood expression of CD59 - an integrative analysis PLoS Comput. Biol. 13 2017 e1005766
53 Taylor D.L. Jackson A.U. Narisu N. Hemani G. Erdos M.R. Chines P.S. Swift A. Idol J. Didion J.P. Welch R.P. Integrative analysis of gene expression, DNA methylation, physiological traits, and genetic variation in human skeletal muscle Proc. Natl. Acad. Sci. USA 116 2019 10883 10888 31076557
54 Theusch E. Chen Y.-D.I. Rotter J.I. Krauss R.M. Medina M.W. Genetic variants modulate gene expression statin response in human lymphoblastoid cell lines BMC Genom. 21 2020 555
55 Peng S. Deyssenroth M.A. Di Narzo A.F. Cheng H. Zhang Z. Lambertini L. Ruusalepp A. Kovacic J.C. Bjorkegren J.L.M. Marsit C.J. Genetic regulation of the placental transcriptome underlies birth weight and risk of childhood obesity PLoS Genet. 14 2018 e1007799
56 Pashos E.E. Park Y. Wang X. Raghavan A. Yang W. Abbey D. Peters D.T. Arbelaez J. Hernandez M. Kuperwasser N. Large, Diverse Population Cohorts of hiPSCs and Derived Hepatocyte-like Cells Reveal Functional Genetic Variation at Blood Lipid-Associated Loci Cell Stem Cell 20 2017 558 570.e10 28388432
57 Panopoulos A.D. D'Antonio M. Benaglio P. Williams R. Hashem S.I. Schuldt B.M. DeBoever C. Arias A.D. Garcia M. Nelson B.C. iPSCORE: A Resource of 222 iPSC Lines Enabling Functional Characterization of Genetic Variation across a Variety of Cell Types Stem Cell Rep. 8 2017 1086 1100
58 Hoffman G.E. Bendl J. Voloudakis G. Montgomery K.S. Sloofman L. Wang Y.C. Shah H.R. Hauberg M.E. Johnson J.S. Girdhar K. CommonMind Consortium provides transcriptomic and epigenomic data for Schizophrenia and Bipolar Disorder Sci. Data 6 2019 180 31551426
59 Guelfi S. D'Sa K. Botía J.A. Vandrovcova J. Reynolds R.H. Zhang D. Trabzuni D. Collado-Torres L. Thomason A. Quijada Leyton P. Regulatory sites for splicing in human basal ganglia are enriched for disease-relevant information Nat. Commun. 11 2020 1041 1116 32098967
60 Steinberg J. Southam L. Roumeliotis T.I. Clark M.J. Jayasuriya R.L. Swift D. Shah K.M. Butterfield N.C. Brooks R.A. McCaskie A.W. A molecular quantitative trait locus map for osteoarthritis Nat. Commun. 12 2021 1309 33637762
61 Young A.M.H. Kumasaka N. Calvert F. Hammond T.R. Knights A. Panousis N. Park J.S. Schwartzentruber J. Liu J. Kundu K. A map of transcriptional heterogeneity and regulatory variation in human microglia Nat. Genet. 53 2021 861 868 34083789
62 Bossini-Castillo L. Glinos D.A. Kunowska N. Golda G. Lamikanra A.A. Spitzer M. Soskic B. Cano-Gamez E. Smyth D.J. Cattermole C. Immune disease variants modulate gene expression in regulatory CD4+ T cells Cell Genom. 2 2022 None
63 Momozawa Y. Dmitrieva J. Théâtre E. Deffontaine V. Rahmouni S. Charloteaux B. Crins F. Docampo E. Elansary M. Gori A.S. IBD risk loci are enriched in multigenic regulatory modules encompassing putative causative genes Nat. Commun. 9 2018 2427 29930244
64 Fairfax B.P. Makino S. Radhakrishnan J. Plant K. Leslie S. Dilthey A. Ellis P. Langford C. Vannberg F.O. Knight J.C. Genetics of gene expression in primary immune cells identifies cell type-specific master regulators and roles of HLA alleles Nat. Genet. 44 2012 502 510 22446964
65 Fairfax B.P. Humburg P. Makino S. Naranbhai V. Wong D. Lau E. Jostins L. Plant K. Andrews R. McGee C. Knight J.C. Innate immune activity conditions the effect of regulatory variants upon monocyte gene expression Science 343 2014 1246949
66 Kasela S. Kisand K. Tserel L. Kaleviste E. Remm A. Fischer K. Esko T. Westra H.J. Fairfax B.P. Makino S. Pathogenic implications for autoimmune mechanisms derived by comparative eQTL analysis of CD4+ versus CD8+ T cells PLoS Genet. 13 2017 e1006643
67 Gilchrist J.J. Makino S. Naranbhai V. Sharma P.K. Koturan S. Tong O. Taylor C.A. Watson R.A. de Los Aires A.V. Cooper R. Natural Killer cells demonstrate distinct eQTL and transcriptome-wide disease associations, highlighting their role in autoimmunity Nat. Commun. 13 2022 4073 35835762
68 Di Angelantonio E. Thompson S.G. Kaptoge S. Moore C. Walker M. Armitage J. Ouwehand W.H. Roberts D.J. Danesh J. INTERVAL Trial Group Efficiency and safety of varying the frequency of whole blood donation (INTERVAL): a randomised trial of 45 000 donors Lancet 390 2017 2360 2371 28941948
69 Zhao H. Sun Z. Wang J. Huang H. Kocher J.P. Wang L. CrossMap: a versatile tool for coordinate conversion between genome assemblies Bioinformatics 30 2014 1006 1007 24351709
70 Kolberg L. Raudvere U. Kuzmin I. Adler P. Vilo J. Peterson H. g:Profiler-interoperable web service for functional enrichment analysis and gene identifier mapping (2023 update) Nucleic Acids Res. 51 2023 W207 W212 37144459
71 Murphy A.E. Schilder B.M. Skene N.G. MungeSumstats: a Bioconductor package for the standardization and quality control of many GWAS summary statistics Bioinformatics 37 2021 4593 4596 34601555
72 Chick J.M. Munger S.C. Simecek P. Huttlin E.L. Choi K. Gatti D.M. Raghupathy N. Svenson K.L. Churchill G.A. Gygi S.P. Defining the consequences of genetic variation on a proteome-wide scale Nature 534 2016 500 505 27309819
73 Gonçalves E. Fragoulis A. Garcia-Alonso L. Cramer T. Saez-Rodriguez J. Beltrao P. Widespread Post-transcriptional Attenuation of Genomic Copy-Number Variation in Cancer Cell Syst. 5 2017 386 398.e4 29032074
