Most autoimmune disease risk effects identified by genome-wide association studies (GWAS) localize to open chromatin with gene regulatory activity. GWAS loci are also enriched for expression quantitative trait loci (eQTLs), suggesting that most risk variants alter gene expression 1,2 . However, because causal variants are difficult to identify and cis-eQTLs occur frequently, it remains challenging to identify specific instances of disease-relevant changes to gene regulation. Here, we use a novel joint likelihood framework with higher resolution than previous methods to identify loci where autoimmune disease risk and an eQTL are driven by a single, shared genetic effect. Using eQTLs from three major immune subpopulations, we find shared effects in only ~25% of loci. Thus, we uncover a fraction of gene regulatory changes as strong mechanistic hypotheses for disease risk, but conclude that most risk mechanisms likely do not involve changes to basal gene expression. The autoimmune and inflammatory diseases (AID) are heritable, complex diseases where loss of tolerance to self-antigens results in either systemic or tissue-specific immune attack 3,4 . GWAS have identified hundreds of genomic regions mediating risk to several AID. These associations are primarily non-coding: lead GWAS SNPs are more likely to be associated with expression levels of neighboring genes than expected by chance 12,13 , and the same lead SNPs are enriched in regulatory regions marked by chromatin accessibility and modification 1,14 . Fine-mapping reveals enrichment of AID-associated variants in enhancer elements active in stimulated T cell subpopulations 15 , with heritability strongly enriched in such regulatory regions 16,17 . Collectively, these strands of evidence suggest that the majority of disease risk is mediated by changes to gene regulation in specific cell subpopulations. However, these bulk analyses do not formally assess whether expression levels and disease risk can be attributed to a single underlying variant or to independent effects in a locus 18,19 . Though several methods have been developed to assess these alternatives using eQTL data 20–23 , they show limited resolution to detect cases where distinct disease and eQTL causal variants are in linkage disequilibrium. Here, we present an approach to test if a GWAS risk association and an eQTL are driven by the same underlying genetic effect, accounting for the LD between causal variants. Using data from ImmunoChip studies of seven AID comprising >180,000 samples in total (Supplementary Table 1), we test if associations in 272 known risk loci are consistent with cis-eQTL for genes in each region, measured in three relevant immune cell populations: lymphoblastoid cell lines (LCLs), CD4+ T cells and CD14+ monocytes 24,25 . When associations to two traits – here, disease trait and eQTL – are driven by the same underlying causal variant, the joint evidence of association should be maximized at the markers in tightest LD with the causal variant 19,26 . Here, we directly evaluate this joint likelihood (Supplementary Figure 1), unlike previous approaches that look for similarities in the shape of the association curve over multiple markers 20,21,27,28 . When the underlying causal effect is shared, joint likelihood is maximized when we model the same causal variant in both traits; conversely, when the underlying causal variants are different, we expect maximum joint likelihood when we model their closest proxies. We empirically derive the null distribution of the joint likelihood ratio statistic by comparing disease associations to permuted eQTL data(see Methods, Supplementary Figure 2 and Supplementary Notes). We thus directly evaluate whether two associations in the same locus, observed in different cohorts, are due to the same underlying effect. To assess the performance of our method, we benchmarked it against three recently reported methods: coloc 20 , a well-calibrated Bayesian framework that considers spatial similarities in association data across sets of markers; gwas-pw 29 , which extends this idea to hierarchical priors and optimizes model parameters; and HEIDI/SMR 22 , which applies Mendelian randomization between traits. We simulated pairs of case-control cohorts with either the same or distinct causal variants driving association, and find that our approach shows the best overall performance (Supplementary Tables 2 and 3). When independent causal variants (i.e. not in LD) drive GWAS and eQTL associations, our own method, coloc and gwas-pw all had excellent performance. As the LD between the causal variants increases, our method shows the best performance, maintaining high resolution even when the underlying causal variants are in strong LD (AUC = 0.883 when 0.7 0.8), where the false positive rate increases dramatically. We also have limited resolution when the eQTL effect is very weak (p > 0.01, Supplementary Figures 12–15). Thus, within these limits, we can accurately detect cases of shared genetic effects between two traits. To dissect AID risk loci, we first identified densely genotyped ImmunoChip loci showing genome-wide significant association, excluding the Major Histocompatibility Locus due to the extensive LD structure in the region (immunobase.org; Table 1). We next identified genes in a 1Mb window centered on the most associated variant in each locus. Consistent with previous observations that eQTLs are frequently found in GWAS loci, we found that 260/272 loci had at least one gene with an eQTL (p 75% of tested disease-eQTL pairs appear associated to distinct genetic variants in the same locus (Figure 1). We sought to explain this lack of overlap between disease associations and eQTLs, despite their frequent co-occurrence in the same loci. In particular, although our method showed good performance in simulated data (Supplementary Figure 4), we remained concerned that this lack of overlap may be due to low statistical power in the eQTL data, which come from cohorts of limited sample size. However, we find that even amongst the most strongly supported eQTLs (nominal p 10% in CEU. We generate genotypes for the eQTL cohorts using HapGen2 with a null effect size, then simulate a quantitative phenotype using the allelic mean difference model implemented in GCTA 36 with effect sizes of 0.05, 0.1, or 0.2 in cis-heritability (h2 ). In each locus, we generated five replicate eQTL cohorts for H0 and H1 each; for H2, we generated a single cohort up to five distinct causal variants per locus. We used plink to calculate the genetic association with disease and expression phenotypes in logistic and linear regression models, respectively, after filtering out SNPs with MAF 10−5 for disease cohorts and association p > 0.01 for eQTL cohorts). In addition, as expected from the coalescent forward simulation model on which HapGen2 is based, a fraction of our simulated cohorts showed maximal association to a SNP in low LD with the causal variant we had specified (r2 5%). We removed pseudogenes and transcripts without assigned gene symbols from the expression data, and calculated association statistics by linear regression of genotype on expression levels, including three population principal components to control for structure 37,38 . For CD4+ and CD14+, we regressed normalized expression levels for European Americans (n=213 and 211, respectively) on similarly QCed imputed allele dosages. For all cell types, we generated adaptive permutation statistics from 103 up to 106 iterations, using all covariates 37 . Joint likelihood mapping (JLIM) To test the hypothesis that association signals for two traits are driven by the same causal variant, we contrasted the joint likelihood of observed association statistics under the assumption of same compared to distinct causal variant. Due to limited genetic resolution, distinct causal variants were defined by separation in LD space by r 2 θ . We derived N θ 1 from the reference LD panel and N θ 2 directly from the genotypes of secondary trait cohort. We used disease outcome as primary trait, leveraging the larger sample size and dense genotyping, and gene expression as secondary trait, taking advantage of the availability of individual genotype data. The likelihood of causal association was calculated by approximating the local LD structure with pairwise correlation similarly as Kichaev et al. and Hormozdiari et al. Briefly, when SNP c is the only causal variant in the locus with non-centrality λc , association static z i of non-causal SNP i follows a normal distribution N(ri,c λc, 1), where ri,c is LD between SNPs i and c measured in pairwise Pearson correlation of genotypes. In general, when association statistics Z = (z 1, z 2,…zM ) T are provided for all M SNPs in the analysis window, the likelihood of SNP i being the causal variant with non-centrality λi is: L ( Z ; λ i ≠ 0 ) = ϕ MVN ( Z ; mean = ∑ ( λ i ∘ C ) , var = ∑ ) where ϕMVN is the multivariate normal density function, C is an incident vector with Ck = 1 if and only if K = I, Σ is a M × M local LD matrix defined by pairwise Pearson correlation between genotypes, and ∘ is element-wise multiplication 39 . Since we do not know the true non-centrality of causal variant, we estimated the profile likelihood, which simplifies to a closed form 40 : log L ( Z ; λ i mle ) = 1 2 ( - Z T ∑ - 1 Z + z i 2 ) - 1 2 log ( ( 2 π ) M ∣ ∑ ∣ ) with λ i mle = z i . Thus, given association statistics for primary and secondary traits, Z = (z 1, z 2,…zM ) T and W = (w 1, w 2,…wM ) T , the test statistic Λ simplifies to: Λ = ∑ i ∈ N θ 1 ( m ∗ ) e 1 2 ( z i 2 - z m ∗ 2 ) · ( w i 2 - max j ∉ N θ 2 ( i ) w j 2 ) The p-value of joint likelihood is estimated by permuting phenotypes of secondary traits as under the trivial null hypothesis that that there is no casual variant for secondary trait in the locus (H0). With respect to the more likely null that distinct causal variants underlie association signals of two traits (H2), we can show that asymptotically as the non-centrality of causal variant increases, p-values estimated from H0 behave conservatively with respect to H2 (Supplementary Notes): P JLIM = P ( Λ ≥ l ∣ H 0 ) ≥ P ( Λ ≥ l ∣ H 2 ) Thus, with large enough sample or effect sizes, joint likelihood test against H0 will also reject H2 in favor of alternative hypothesis of shared causal variant (H1). Further, to evaluate whether this property holds for practical non-centrality values, we examined our negative controls simulating H2, specifically, if PJLIM was highly shifted toward 1.0 (Supplementary Figure 2) and similar or larger than empirically estimated false positive rates as expected (Supplementary Table 3 ≤ 0.05). For both simulated and real GWAS data, we applied JLIM to SNPs with data for both primary and secondary traits, present in the reference LD panels, and within 100kb of the most associated marker to disease (“lead SNP”). In ImmunoChip data, the analysis windows were further confined by the boundaries of the dense genotyping intervals. We compared each lead SNP to eQTL data for all genes with a transcription start sites (TSS) up to 1Mb from the lead SNP, and an eQTL association p 0.1 to both). For the reference LD panel, we used the base haplotypes of HapGen simulation for simulated datasets, and non-Finnish European samples (n=404) of the 1000 Genomes Project (phase 3, release 2013/05/02) for ImmunoChip loci. We corrected for multiple tests using false discovery rate (FDR) levels and Bonferroni correction. The FDR was calculated separately for specific disease and cell type combination as: FDR ( p ) = p N # { P JLIM ≤ p } where p is a JLIM p-value cut-off, and N is the number of all tested disease lead SNP-eQTL candidate gene combinations. The FDR was calculated for each cell type since the distribution of JLIM p-values can vary depending on the disease relevance of cell type. To provide a list of higher confidence hits in each disease, we also applied the Bonferroni correction to nominal JLIM p-values for the number of tests across all three cell types. Benchmark comparison We used our simulations to compare the performance of our method (here abbreviated as JLIM, for joint likelihood mapping) to three existing methods: Bayesian coloc, gwas-pw, and SMR/HEIDI. We ran coloc (version 2.3-1) using default parameter settings with the colocalization prior p12 set to 10−6. We followed the authors’ recommendation to use beta and variance of beta for the case/control cohorts as summary statistics. For quantitative trait cohorts, we also provided in-sample minor allele frequencies. We applied gwas-pw (version 0.21) with default parameters. All simulated disease-eQTL pairs were combined into a single batch and analyzed together so that gwas-pw can optimize the model parameters. We ran SMR/HEIDI (version 0.64) with default parameters except the p-value threshold to select the top associated eQTL (peqtl-smr, default 5 × 10−8), which was relaxed to 0.01 in order to enable the test on simulated disease-eQTL pairs with weak eQTL association. Tests producing significant heterogeneity by HEIDI (pHEIDI < 0.05) were called negative regardless of pSMR values since they are likely to harbor distinct causal variants between disease and eQTL. For coloc and gwas-pw, predictions were made based only on reported posterior probability of colocalization (PP4 and PP3, respectively) although they report posteriors for other competing models. For overall performance comparison, we evaluated the area under the receiver operator curve (ROC; Supplementary Figures 3, p and 4) using H1 as known positives and H0 and H2 as known negatives. For SMR/HEIDI, the sensitivity did not reach 1.0 even at the specificity of 0.0 since the method called significant heterogeneity on 15% of H1 (pHEIDI < 0.05). The sensitivity and specificity were also compared at p-value cut-offs of 0.01 and 0.05 (Supplementary Tables 2 and 3). For coloc and gwas-pw, the posterior probability cut-offs equivalent to p-values were determined from the false positive rates on null simulation of no eQTL. Bayesian coloc on real data As ImmunoChip data is only available as summary statistics, we used ran coloc 20 with the p-values of association the minor allele frequencies from non-Finnish Europeans from the 1000 Genomes Projectfor disease cohorts, and with quantitative beta and variance of beta calculated onfor eQTL association datacohorts. We also provided to coloc the minor allele frequencies of non-Finnish Europeans from the 1,000 Genomes Project 35 . The and a colocalization prior p12 was set to= 10−6 6, and the prediction was made at PP4 ≥ 0.75 for higher confidence (Supplementary Table 4). We did not consider the type 1 diabetes data, where case/control sample size is limited after excluding affected sib pair data. Estimating the number of disease GWAS loci with consistent eQTL effects We expect JLIM p-values to follow a bimodal distribution with modes close to zero and one when the data support a model of shared or distinct causal effects, respectively. Conversely, under the null model of no cis-eQTL association, we expect a uniform p-value distribution. We can thus estimate the proportion of disease-eQTL pairs belonging to the null π 0, same π 1 and distinct π 2 causal variant models from the observed p-value distribution 41 (Supplementary Figures 16–18). To assess if the strength of the eQTL association influences the likelihood of identifying a shared causal variant, we calculate these proportions for subsets of trait pairs defined by minimum eQTL p-value. In each bin, we identified the limits of the uniform portion of the distribution γ 1 and γ 2 and estimate π 0, π 1 and π 2 as: π 0 = # { γ 1 < P JLIM < γ 2 } ( γ 2 - γ 1 ) N π 1 = # { P JLIM ≤ γ 1 } N - γ 1 π 0 π 2 = # { P JLIM ≥ γ 2 } N - ( 1 - γ 2 ) π 0 To estimate the number of disease GWAS loci that can be explained by consistent effect of same causal variant on disease and eQTL (denoted by 𝒞 below), we incrementally relaxed the p-value cut-offs of JLIM and examined the trends of the number of disease loci with at least one JLIM hit and subtracted the expected number of false positive loci (Figure 1 and Supplementary Figure 19). Specifically, at each JLIM p-value cutoff pi , we successively calculated 𝒞(pi ): C ( p i ) = C ( p i - 1 ) + ∣ D ( p i ) ∩ D ( p i - 1 ) c ∣ - ∑ d ∈ D ( p i - 1 ) c E ( d , p i ) + ∑ d ∈ D ( p i - 1 ) c E ( d , p i ) where pi −1 < pi with p 0 = 0, 𝒟(p) is the set of disease GWAS loci with at least one eQTL gene in any cell type passing the JLIM p-value cut-off p, and ε (d,p) is the probability that disease GWAS locus d has a false positive eQTL gene passing the JLIM p-value cutoff p. We estimated the lower and upper bounds of ε (d,p) using the Monte Carlo method by randomly selecting false positive eQTL genes within the locus d at rates of (1 − π 1) · lb or (1 − π 1) · ub over 1,000 iterations. The lb and ub are the lower and upper bounds of false positive rate of JLIM against true null. Note that π 1 and lb depend on the cell type and strength of eQTL association. As the true null is mixture of two nulls, H 0 and H 2, the false positive rate of JLIM against true null P(Λ ≥ l|H 0 ∪ H 2) can be bounded by using the following decomposition: P ( Λ ≥ l ∣ H 0 ∪ H 2 ) = P ( Λ ≥ l ∣ H 0 ) P ( H 0 ) P ( H 0 ) + P ( H 2 ) + P ( Λ ≥ l ∣ H 2 ) P ( H 2 ) P ( H 0 ) + P ( H 2 ) While the false positive rate under distinct null P(Λ ≥ l|H 2) is difficult to estimate, it is non-negative by definition and asymptotically bounded by permutation p-value P(Λ ≥ l|H 0), i.e. PJLIM , as the non-centrality of causal variant increases. Therefore, we took: u b = P JLIM l b = P JLIM π 0 π 0 + π 2 = P JLIM π 0 1 - π 1 and estimated the bounds of locus-level false positive rates ε(d,p) and number of disease loci with consistent effects 𝒞(pi ). Supplementary Material 1
