A fast and efficient colocalization algorithm for identifying shared genetic risk factors across multiple traits. — TEI-tagged view

Christopher N. FoleypersName, James R. StaleypersName, Philip G. BreenpersName, Benjamin B. SunpersName, Paul D. W. KirkpersName, Stephen BurgesspersName, Joanna M. M. HowsonpersName (2021date). “A fast and efficient colocalization algorithm for identifying shared genetic risk factors across multiple traits..” Nature Communicationstitle·j.

Object of study — multi-trait statistical colocalizationterm·mth of genetic association signals

Research question — how to detect colocalizationterm·key across large numbers of traits without enumerating every causal configuration

Contribution — HyPrColocterm·mth, a deterministic Bayesianterm·key algorithm that approximates multi-trait colocalizationterm·key in linear time

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

A fast and efficient colocalizationterm·key algorithm for identifying shared genetic risk factors across multiple traits.

Abstract Genome-wide association studiesterm·mth (GWASterm·mth) have identified thousands of genomic regions affecting complex diseases. The next challenge is to elucidate the causal genes and mechanisms involved. One approach is to use statistical colocalizationterm·mth to assess shared genetic aetiology across multiple related traits (e.g. molecular traits, metabolic pathways and complex diseases) to identify causal pathways, prioritize causal variantsterm·key and evaluate pleiotropyterm·key. We propose HyPrColocterm·mth (Hypothesis Prioritisation for multi-trait Colocalizationterm·mth), an efficient deterministic Bayesianterm·key algorithm using GWAS summary statisticsterm·mth that can detect colocalizationterm·key across vast numbers of traits simultaneously (e.g. 100 traits can be jointly analysed in around 1 s). We perform a genome-wide multi-trait colocalizationterm·key analysis of coronary heart diseaseterm·dom (CHDterm·dom) and fourteen related traits, identifying 43 regions in which CHDterm·dom colocalized with ≥1 trait, including 5 previously unknown CHDterm·dom loci. Across the 43 loci, we further integrate gene and protein expression quantitative trait loci to identify candidate causal genes.

Introduction Genome-wide association studiesterm·mth (GWASterm·mth) have identified thousands of genomic loci associated with complex traits and diseases ( https://www.ebi.ac.uk/gwasterm·mth/ ). However, identification of the causal mechanisms underlying these associations and subsequent biological insights have not been as forthcoming, due to issues such as linkage disequilibriumterm·dom (LDterm·dom) and incomplete genomic coverage. One approach to aid biological insight following GWASterm·mth is to make use of functional data. For example, candidate causal genes can be proposed when the overlap in association signals between a complex trait and functional data (e.g. gene expression) is a consequence of both traits sharing a causal variantterm·key, i.e. the association signals for both traits colocalize. The abundance of significant associations identified by GWASterm·mth means that chance overlap between association signals for different traits is likely 1 . Consequently, overlap does not by itself allow us to identify causal variantsterm·key 1 , 2 . Statistical colocalizationterm·mth methodologies seek to resolve this. By constructing a formal statistical model, colocalizationterm·key approaches have been successful in identifying whether a molecular trait (e.g. gene expression) and a disease trait share a causal variantterm·key in a genomic region 3 – 7 , and potentially prioritise a candidate causal gene. Recently it has been proposed that colocalizationterm·key methodologies can be further enhanced by integrating additional information from multiple intermediate traits linked to disease, e.g. protein expression, metabolite levels 8 . The underlying hypothesis of multi-trait colocalizationterm·key is that if a variant is associated with multiple related traits then this provides stronger evidence that the variant may be causal 8 . Thus, multi-trait colocalizationterm·key aims to increase power to identify causal variantsterm·key. We show that by using multi-level functional datasets in this way can reveal candidate causal genes and pathways underpinning complex disease.

A number of statistical methods have been developed to assess whether association signals across a pair of traits colocalize 3 – 7 . Some methods assess colocalizationterm·key between a pair of traits using individual participant dataterm·mth 9 , 10 , limiting their applicability. In contrast, the COLOC algorithm uses GWAS summary statisticsterm·mth 2 . This approach works by systematically exploring putative causal configurations, where each configuration locates a causal variantterm·key for one or both traits, under the assumption that there is at most one causal variantterm·key per trait. COLOC was recently extended to the multi-trait framework, MOLOC 8 . The authors achieved a 1.5-fold increase in candidate causal gene identification when a third relevant trait was included in colocalizationterm·key analyses relative to results from two traits. However, the approach is computationally impractical beyond 4 traits due to prohibitive computational complexity arising from the exponential growth in the number of causal configurations that must be explored with each additional trait analysed.

Here we present a computationally efficient method, hypothesis prioritisation for multi-trait colocalizationterm·mth (HyPrColocterm·mth), to identify colocalized association signals using summary statisticsterm·key on large numbers of traits. The approach extends the underlying methodology of COLOC and MOLOC. Our major result is that the posterior probability of colocalizationterm·key at a single causal variantterm·key can be accurately approximated by enumerating only a small number of putative causal configurations. Moreover, HyPrColocterm·mth identifies subsets (which we refer to as clusters) of traits which colocalize at distinct causal variantsterm·key in the genomic locus by employing a novel branch and bound divisive clustering algorithm. We show that the multi-trait clustering method of HyPrColocterm·mth has several performance advantages over alternative colocalizationterm·key approaches and apply HyPrColocterm·mth genome-wide to coronary heart diseaseterm·dom (CHDterm·dom) and many related traits 11 , 12 , identifying known and previously unknown candidate CHDterm·dom genetic risk loci with colocalized associations across these traits.

Results Overview HyPrColocterm·mth is a Bayesianterm·key method for identifying shared genetic associations between complex traits in a particular gene region using summary GWASterm·mth results. HyPrColocterm·mth provides two principal novelties: (i) Efficient computation of the posterior probability that all m traits share a causal variantterm·key (which we refer to as the posterior probability of full colocalizationterm·mth, PPFC); and (ii) partitioning of traits into clusters, such that each cluster comprises traits sharing a causal variantterm·key. HyPrColocterm·mth only requires regression coefficients and their corresponding standard errors from summary GWASterm·mth (for binary traits these can be on the log-odds scale, see Methods). The approach makes three key assumptions: (i) for non-independent studies, that the GWASterm·mth results are from the same underlying population, i.e. that the LDterm·dom pattern is the same across studies, (ii) that there is at most one causal variantterm·key in the genomic region for each trait (we assess limitations of this assumption when there are multiple causal variantsterm·key below), and (iii) that these causal variantsterm·key are either directly typed or well imputed in all of the GWASterm·mth datasets 2 , 8 .

Description of the HyPrColocterm·mth method We define a putative causal configuration matrix S to be a binary m × Q matrix, where m is the number of traits and Q is the number of variants. To increase the probability of identifying any underlying causal variantterm·key(s) in the region, the number of SNPs Q included in analyses should be maximised, i.e. the region should be well imputed. S ij is 1 if the j th variant is causal for the i th trait and 0 otherwise ( Supplementary Information ). A hypothesis uniquely identifies traits which share a causal variantterm·key, traits which have distinct causal variantsterm·key and traits which do not have a causal variantterm·key. Except for the null hypothesis ( H 0 ) of no causal variantterm·key for any trait, hypotheses such as H m : all m traits share a causal variantterm·key correspond to multiple configuration matrices, S (Fig. 1 ). By considering the set of configurations to which a hypothesis corresponds, the posterior odds of the hypothesis against the null hypothesis can be computed. For example, let S m denote the set of configurations for hypothesis H m and S 0 denote the single configuration for H 0 , then the posterior odds for the hypothesis that all traits colocalize to a single causal variantterm·key is given by, 1 P ( H m ∣ D ) P ( H 0 ∣ D ) = ∑ S ∈ S m P ( D ∣ S ) P ( D ∣ S 0 ) × p S p S 0 Fig. 1 Colocalizationterm·key hypotheses and causal configurations. Statistical colocalizationterm·mth hypotheses and examples of their associated SNP configurations that allow for at most one causal variantterm·key for each of m traits in a region containing Q genetic variants. For clarity, the hypotheses and a single configuration associated with each hypothesis are shown for m ≥ 4 traits, but the column totals Bell ( m + 1) and ( Q + 1) m are correct for m ≥ 2.

where D represents the combined trait data, the first term in the summation is a Bayes factor and the second term is a prior odds 2 , 8 . To identify a candidate causal variantterm·key across the m traits, i.e. to perform multi-trait fine-mappingterm·key, we locate the configuration S * satisfying max S ∈ S m P S D = P S * D . If the summary data for the genetic associations between traits are independent, then the Bayes factor for each configuration S can be computed by combining Wakefield’s approximate Bayes factors 13 for each trait in the configuration (‘Methods’). If the summary data between traits are correlated because a subset of the participant data was used in at least two of the GWASterm·mth analyses, then an extension to Wakefield’s approximate Bayes factors, which jointly models the trait associations, can be employed (‘Methods’). For a given hypothesis H and set of corresponding configurations S H , the prior probability of configuration S , p ( S ), can either be equal for all S ∈ S H , or can be defined as a product of variant-level priors (‘Methods’). Our variant-level prior extends that of COLOC 2 and MOLOC 8 to a framework that is suitable for the analysis of large numbers of traits. We adopt an approach which requires the specification of a partition of the traits into clusters, together with two interpretable parameters: p , the probability that a variant is causal for one trait; and p c , the conditional probability that a variant is causal for a second trait given it is causal for one trait (‘Methods’). As it will be helpful later, we refer to p c as the conditional colocalizationterm·key prior. COLOC 2 requires specification of three prior parameters { p 1 , p 2 , p 12 } and, while the scope of the configuration priors in HyPrColocterm·mth is different for more than a pair of traits, it is instructive to note that p ≡ pi , for i ∈ { 1 , 2 } , and p c ≡ p 12 p 12 + p 1 when m = 2. To help users of the COLOC 2 software, our software allows users to specify the parameter p and one of either (i) p c ; or (ii) p 12 , from which p c is computed. For simplicity and as a conservative measure, we assume a priori that the genetic association probability p and the conditional colocalizationterm·key probability p c are equal for all traits. This approach allows sensitivity analyses assessing robustness of posterior inference to be routinely performed. However, it implicitly assumes traits are a priori exchangeable, e.g. assumes p 1 = p 2 ; this is supported across a range of designs (case/control or quantitative trait) but may lead to poorer performance in specific datasets 14 .

Efficient computation of the posterior probability of full colocalizationterm·mth (PPFC) For a pre-specified genomic region comprising Q variants, the aim is to evaluate the PPFC , P ( H m | D ), that all m traits share a causal variantterm·key within that region, given the summarized data D . According to Bayes’ rule, this is given by: 2 P P F C : P H m ∣ D = ∑ S ∈ S m P ( D ∣ S ) × p S p ( D ) .

Brute-force computation of the denominator, p ( D ), requires the exhaustive enumeration of ( Q + 1) m causal configurations, which is computationally prohibitive for m > 4, e.g. MOLOC 8 . HyPrColocterm·mth overcomes this challenge by approximating p ( D ) in a way that is both computationally efficient, i.e. has fast computational time, and tightly bounds the approximation error.

As we show in the ‘Methods’, the PPFC can be approximated as 3 P P F C ^ = P R P A , where P R , P A > 0 are rapidly computable values that quantify the probability that two criteria necessary for colocalizationterm·key are satisfied (Fig. 2 ). The first of these criteria is that all the traits must share an association with one or more variants within the region. P R , which we refer to as the regional association probability, is the probability that this criterion is satisfied. By itself, this criterion does not guarantee that there is a single causal variantterm·key shared by all traits, because it could be the case that two or more traits have distinct causal variantsterm·key in strong LDterm·dom with one another. To safeguard against this, we have a second criterion that ensures the shared associations between all traits are owing to a single shared putative causal variantterm·key. P A is the probability that this second criterion is satisfied. We refer to P A as the alignment probabilityterm·mth as it quantifies the probability of alignment at a single causal variantterm·key between the shared associations. Both P R and P A have linear computational cost in the number of traits m , making a calculation of P P F C ^ possible when analysing vast numbers of traits. If the first criterion is satisfied, but the second is not, this may be because it is possible to partition the traits into clusters, such that each cluster has a distinct causal variantterm·key. HyPrColocterm·mth additionally seeks to identify these clusters. Fig. 2 Illustration of the HyPrColocterm·mth approximation. We illustrate the HyPrColocterm·mth approach with m = 2 traits. Statistical colocalizationterm·mth between traits which do not share an association region, i.e. do not have shared genetic predictors, is not possible (no colocalizationterm·key criteria satisfied). However, traits which do (satisfying criterion 1) possess the possibility. HyPrColocterm·mth first assesses evidence supporting all m traits sharing an association region, which quickly identifies utility in a colocalizationterm·key framework. HyPrColocterm·mth then assesses whether any shared association region is due to colocalizationterm·key between the traits (criteria 1 and 2) or due to a region of strong LDterm·dom between two distinct causal variantsterm·key, one for each trait (criterion 1 only). Results from these two calculations are combined to accurately approximate the PPFC .

Identification of clusters of colocalized traits If P P F C ^ falls below a threshold value, τ , we reject the hypothesis H m that all m traits colocalize to a shared causal variantterm·key. In practice, this threshold is specified by defining separate thresholds, P R * and P A * , for P R and P A , such that τ = P R * P A * (‘Methods’). If H m is rejected, HyPrColocterm·mth seeks to determine if there are values ℓ < m such that H ℓ cannot be rejected; i.e. if there exist subsets of the traits such that all traits within the same subset colocalize to a shared causal variantterm·key. Starting with a single cluster containing all m traits, our branch and bound divisive clustering algorithm (Supplementary Figs. S1a, b ) iteratively partitions the traits into larger numbers of clusters, stopping the process of partitioning a cluster of two or more traits when all traits in a cluster satisfy both P R > P R * and P A > P A * . The process of partitioning a cluster into two smaller clusters is performed using one of two criteria: (i) regional ( P R ) or (ii) alignment ( P A ) selection (‘Methods’ and Supplementary Note ). For k ≤ m traits in a cluster, the regional selection criterion has O k Q computational cost and is computed from a collection of hypotheses that assume not all traits in a cluster colocalize because one of the traits does not have a causal variantterm·key in the region. The alignment selection criterion has O k Q 2 computational cost and is computed from hypotheses that assume not all traits in a cluster colocalize because one of the traits has a causal variantterm·key elsewhere in the region ( Supplementary Note ). By default, the HyPrColocterm·mth software uses the more computationally efficient regional selection criterion to partition a cluster.

Model validation using simulations We created simulated datasets by resampling phased haplotypes from the European samples in 1000 Genomes 15 and for each dataset we randomly selected one of the first 50 regions confirmed to be associated with CHDterm·dom 16 (‘Methods’). For each simulation scenario, 1000 replicates were performed.

Computational efficiency The posterior probability of colocalizationterm·key, across m traits and in a region of Q variants, can be accurately approximated by computing O m Q 2 causal configurations. Figure 3 illustrates this for varying numbers of independent studies and variants, demonstrating a close linear relationship between computation time and the number of traits. Consequently, HyPrColocterm·mth is able to assess 100 traits, in a region of 1000 SNPs, in under 1 second compared to MOLOC which takes approximately 1 hour to analyse five traits. For m ≤ 4, traits the median absolute relative difference between the HyPrColocterm·mth and MOLOC 8 posterior probabilities was found to be ⪅ 0.5% (Fig. 3 ). Fig. 3 Comparison of HyPrColocterm·mth and MOLOC computation time and posterior probability of colocalizationterm·key. (Left panel) Computation time (seconds) for HyPrColocterm·mth (yellow) and MOLOC (blue) to assess full colocalizationterm·key across M ≤ 1000 traits in a region containing Q = 1000 SNPs. MOLOC was restricted to M ≤ 5 traits owing to the computational and memory burden of the MOLOC algorithm when M > 5. When M = 5, we summarise the computation time of MOLOC from 10 datasets - as it took around 1 hour to analyse a single dataset, in all other scenarios performance was summarised from 1000 datasets. Three reference lines are plotted: (i) Bell ( M + 1), which denotes the theoretical cost of exhaustively enumerating all hypotheses; (ii) M 2 , denoting quadratic cost and; (ii) M 1 , denoting the linear complexity of the HyPrColocterm·mth algorithm. (Right panel) Distribution of the posterior probability of colocalizationterm·key between all traits, i.e. the posterior probability of full colocalizationterm·mth (PPFC), using HyPrColocterm·mth (yellow) and MOLOC (blue) across M ∈ {2,3,4} traits. Error bars denote the 1st and 9th deciles and a point denotes the median value. Despite differences in the prior set-up between the methods, the median absolute relative difference between the two posterior probabilities was ≲ 0.005.

Performance of HyPrColocterm·mth to detect multi-trait colocalizationterm·key We used simulated datasets in which all traits colocalize to assess the accuracyterm·mth of HyPrColocterm·mth in detecting colocalizationterm·key across varying numbers of traits and study sample sizes. We simulated independent datasets with sample sizes of 5,000, 10,000, and 20,000 individuals for up to 100 quantitative traits and for which all traits share a single causal variantterm·key explaining either 0.5%, 1% or 2% of trait variance. For each simulated dataset, we used HyPrColocterm·mth to approximate the PPFC. The distribution of PPFC across the simulated datasets was narrower in the analysis of two traits relative to a larger number of traits, as the probability of random misalignment of the lead variant between traits increases as the number of traits increases (top Fig. 4 ). However, the estimated PPFC is always close to 1 for 5, 10 and 20 traits illustrating that the distribution of the estimate is stable across a broad number of traits and sample sizes. For 100 traits there is a small decrease in power due to the growth in the number of hypotheses in which only a subset of the traits colocalize. This is expected when sample size is fixed and the shared causal variantterm·key explains only a small fraction of trait variation for each trait, as combined evidence supporting hypotheses in which a subset of the traits colocalize are eventually greater than evidence supporting full colocalizationterm·key. Fig. 4 Assessment of the HyPrColocterm·mth posterior probability. Simulation results for a sample size N ∈ {5,000, 10,000, 20,000} and a causal variantterm·key explaining {0.5%, 1%, 2%} of variation across m ∈ {2, 5, 10, 20, 100} traits. Presented is the distribution of the HyPrColocterm·mth posterior probability of full colocalizationterm·mth (PPFC) for variant-level priors only (top); the probability of correctly identifying the causal variantterm·key (middle) and; linkage disequilibriumterm·dom between an incorrectly identified causal variantterm·key and the true causal variantterm·key (bottom). Error bars denote the 1st and 9th deciles and a point denotes the median value and performance was summarised from 1000 simulated datasets. Comparing performance across increasing study sample size and variance explained by the causal variantterm·key, power to detect all colocalized traits is reduced when including studies with smaller sample sizesterm·dom (top row), however, including these studies can still boost the probability of correctly identifying the shared causal variantterm·key irrespective of variance explained (middle row).

When at least one trait did not have a causal variantterm·key in the region the false detection rate was negligible. For example, we generated 100 quantitative traits, each from a study with sample size 10,000, in which 99 traits share a causal variantterm·key and the remaining trait had either: (i) a distinct causal variantterm·key or (ii) no causal variantterm·key in the region. In each scenario a causal variantterm·key explained 1% of trait variation. The 1st, 5th (median) and 9th deciles of the PPFC were (4 × 10 −24 , 1 × 10 −17 , 5 × 10 −8 ) in scenario (i) and (0.02, 0.05, 0.10) in scenario (ii). There is a considerable difference between the results from each scenario, but the PPFC is below the threshold for declaring colocalizationterm·key in both situations.

Fine mappingterm·key the causal variantterm·key with HyPrColocterm·mth The proportion which HyPrColocterm·mth correctly identified the true causal variantterm·key increased as the number of colocalized traits included in the analyses increased up to 2-fold, irrespective of sample size and variance explained by the causal variantterm·key (middle Fig. 4 ), highlighting a major benefit of performing multi-trait fine-mappingterm·key. If HyPrColocterm·mth identified a variant that was not the true causal variantterm·key, we computed the LDterm·dom between the true causal variantterm·key and the identified variant. In cases where the identified variant was not the causal variantterm·key, the variant was typically in very strong LDterm·dom (median r 2 ≥ 0.99) with the true causal variantterm·key and for large numbers of traits, i.e. m ≥ 20, with sample size 20,000, the two variants were in perfect LDterm·dom, i.e. r 2 = 1 (bottom Fig. 4 ).

Branch and bound divisive clustering algorithm Here we assess the performance of the branch and bound (BB) divisive clustering algorithm to identify clusters of colocalized traits over a range of scenarios, several specifications of the conditional colocalizationterm·key prior p c and using three classification criteria. The criteria were: accuracyterm·mth, which is an overall measure of the classification of traits into clusters; true positive rateterm·mth (TPRterm·mth) and false positive rateterm·mth (FPRterm·mth), see Methods for more details. We simulated 10 traits from non-overlapping datasets under three scenarios: (i) a single cluster of 10 colocalized traits; (ii) 2 clusters of 3 colocalized traits, the remaining 4 traits do not have a causal variantterm·key (reflecting hypothesis free colocalizationterm·key searches) and (iii) 4 clusters of colocalized traits, comprising 2 clusters of 3 traits and 2 clusters of 2 traits. Scenarios (ii) and (iii) are designed to simultaneously investigate potential false and true positive findings. Each cluster of colocalized traits share a single causal variantterm·key and causal variantsterm·key between clusters are distinct, but can be in perfect LDterm·dom, i.e. r 2 = 1, with one another—we assess results when the single causal variantterm·key assumption is violated later. To mirror real scenarios in which data are taken from studies with different sample sizes, we take the number of individuals in each study ( N i ) as a random draw from the set N i ∈ {1k, 5k, 10k, 15k, 20k}. For comparison, we additionally present results when all studies have a large sample size by also performing an analysis in which N i = 15k for all traits. In all scenarios, the causal variantterm·key for each trait explained 1% of trait variance and the probability parameters were set to P R * = P A * = 0.5 (‘Methods’). Following the approach of Wallace 14 , we assess sensitivity to the choice of colocalizationterm·key prior p c , i.e. (1 − γ ). Across a wide range of simulated data, Wallace 14 demonstrated that setting p 12 = 5 × 10 −6 in COLOC (approximately p c = 0.05 in HyPrColocterm·mth) was generally a robust choice. Starting from this value, we evaluated results with more conservative choices of p c by performing three separate analyses for each dataset using p c ∈ {0.05, 0.02, 0.01}, equivalent to p 12 ≈ {5 × 10 −6 , 2 × 10 −6 , 1 × 10 −6 } with p = 10 −4 fixed 14 , in order to identify a robust choice of p c . These values can result in substantial differences in the prior probability of colocalizationterm·key as the number of traits in a cluster increases (‘Methods’). For comparison, we compare HyPrColocterm·mth against the alternative of performing pairwise colocalizationterm·key analyses using COLOC 2 , which restricts clusters’ sizes to two traits only. Results are presented in Figs. 5 , 6 and Supplementary Fig. S2 . Fig. 5 Number of clusters of colocalized traits and number of traits within a cluster. Results from the single causal variantterm·key simulation study (c.f. Supplementary Fig. S2 ), presenting a the number of clusters of colocalized traits; and b the number of traits within each cluster identified by HyPrColocterm·mth. Error bars denote the 1st and 9th deciles and a point denotes the median value. Fig. 6 Performance of the BB clustering algorithm when excluding clusters of colocalized traits with lower posterior probability. In each of the three scenarios presented, m = 10 traits with non-overlapping samples were generated, trait sample sizes were drawn randomly from the set N = {1,000, 5,000, 10,000, 15,000, 20,000} and variant-level causal configuration priors were used with three choices of the colocalizationterm·key prior p c ∈ {0.05, 0.02, 0.01}. In scenario (i) there is one cluster of 10 colocalized traits; in scenario (ii) there are 2 clusters of colocalized traits, each comprising of 3 traits, the remaining 4 traits do not have causal variantsterm·key and; in scenario (iii) there are 4 clusters of colocalized traits, 2 clusters of 3 traits and 2 clusters of 2 traits sharing a causal variantterm·key. Traits within a cluster share a single causal variantterm·key and causal variantsterm·key between clusters are distinct, however, a distinct variant can be in perfect LDterm·dom, i.e. r 2 = 1, with another distinct variant. In all scenarios, we present results that passed the posterior probability of colocalizationterm·key P R P A ≥ 0.7. Presented are the classification measures: a accuracyterm·mth; b true positive rateterm·mth; and c the false positive rateterm·mth. See ‘Methods’ for a description of how we define these in the context of clusters of colocalized traits. In d we present the LDterm·dom between the identified causal variantterm·key for each cluster of colocalized traits and the true causal variantterm·key for each cluster. Error bars denote the 1st and 9th deciles and a point denotes the median value. The results highlight that on increasing the posterior threshold from 0.5 (c.f. Supplementary Fig. S2 ) to 0.7, HyPrColocterm·mth’s ability to cluster multiple traits together demonstrably improves accuracyterm·mth and the true positive rateterm·mth relative to pairwise analyses.

We observed that both HyPrColocterm·mth and pairwise COLOC perform reasonably well across all three scenarios. The median accuracyterm·mth and TPRterm·mth is generally ≥0.75, for all three choices of p c , improving to around 1 when the sample size of each study is large (Supplementary Fig. S2a, b ; Supplementary Table S6 )—indicating that including studies with smaller sample sizesterm·dom decreases the TPRterm·mth. Accuracyterm·mth was more sensitive to the choice of p c when all traits colocalized into a single cluster, i.e. scenario (i), relative to scenarios (ii) and (iii) where we observed little sensitivity to p c (Supplementary Fig. S2a ). We noted increased variability in the TPRterm·mth when traits that do not have a causal variantterm·key were included in analyses, i.e. scenario (ii), particularly using the more stringent colocalizationterm·key prior p c = 0.01 (Supplementary Fig. S2b ). The FPRterm·mth was generally low across all scenarios and prior choices: the 1st decile and median values were all zero. However, in scenario (iii), when there are 4 clusters of traits and 4 causal variantsterm·key in the region, the 9 th decile of the FPRterm·mth increased for both methods, from around zero in scenario (ii) up to 0.16, 0.1 and 0.08 when p c was 0.05, 0.02 and 0.01, respectively (Supplementary Fig. S2c ). The increase in FPRterm·mth in scenario (iii) was a consequence of HyPrColocterm·mth occasionally wrongly including an extra trait in one of the clusters (Fig. 5b ), and the pairwise approach overestimating the number of clusters (Supplementary Fig. S1c ). This was because the causal variantsterm·key from distinct clusters were in strong LDterm·dom, i.e. r 2 > 0.95, the FPRterm·mth of both methods reduced when excluding causal variantsterm·key in strong LDterm·dom (Supplementary Fig. S3 ). Over all scenarios, HyPrColocterm·mth regularly identified both the correct number of clusters of colocalized traits in the data (Fig. 5a ) as well as the correct number of colocalized traits within each cluster (Fig. 5b ). The pairwise approach resulted in more variation in the number of clusters identified (Supplementary Fig. S1c ). HyPrColocterm·mth can assign more than a pair of traits to a cluster, allowing information about the location of any shared causal variantterm·key to be borrowed across multiple traits, and therefore performed better at identifying the shared causal variantterm·key (Supplementary Fig. S2d ). HyPrColocterm·mth significantly outperformed the pairwise approach when summarising results from the clusters of colocalized traits whose posterior probability satisfied P R P A > 0.7 (Fig. 6 ; Supplementary Table S7 ). This procedure reflects common practice, as colocalizationterm·key results are generally only reported when the posterior probability of colocalizationterm·key is greater than a threshold value, which we take here to be 0.7. Across all three scenarios, clusters of colocalized traits identified by HyPrColocterm·mth had a median accuracyterm·mth and TPRterm·mth of 1, with little sensitivity to the different choices of colocalizationterm·key prior p c . The FPRterm·mth reduced also, for example in scenario (iii) when p c = 0.01, the 1st, median and 9th deciles of the FPRterm·mth were all zero. The FPRterm·mth reduced for the pairwise approach after thresholding, but the TPRterm·mth reduced as well. In pairwise approaches, a cluster of 3 or more colocalized traits is identified if and only if all pairs of traits colocalize (ideally at the same shared causal variantterm·key), the TPRterm·mth of the pairwise method reduced after thresholding as only some of the pairs of traits passed the posterior threshold which increased the false negative rate. This is a drawback of methods which do not perform multi-trait colocalizationterm·key. We repeated this simulation procedure for 20 traits and the results were similar (Supplementary Fig. S3B ), highlighting the scalability of HyPrColocterm·mth to identify larger clusters of colocalized traits. Overall, across the range of scenarios considered the selection algorithm performed well in terms of sensitivity, specificity and accuracyterm·mth. In many situations there will not be a strong prior belief in a single value for p c . Based on our results and previous investigations 14 , we recommend users set p c = 0.02 and report results from the clusters of colocalized traits which satisfy P R P A > 0.7. Setting p c = 0.02 increased the number of datasets in which clusters satisfying P R P A > 0.7 were identified (‘Methods’) while maintaining a low FPRterm·mth throughout. The HyPrColocterm·mth default p c = 0.02 is equivalent to setting p 12 ≈ 2 × 10 −6 which, for a pair of traits, is slightly more conservative than the recommended value of p 12 = 5 × 10 −6 by Wallace 14 . For more than a pair of traits, however, it can be much more conservative, e.g. setting p c = 0.05 (i.e. p 12 ≈ 5 × 10 −6 ) returns a prior probability of colocalizationterm·key across 10 traits that is around 2000 times larger than when setting p c = 0.02 (i.e. p 12 ≈ 2 × 10 −6 ).

In scenarios (i), (ii) and (iii), HyPrColocterm·mth identified the clusters of colocalized traits on average 50, 30 and 25 times faster than the pairwise COLOC approach, indicating some sensitivity in computational performance to the type of colocalizationterm·key structure present in the data. These figures improved to 200, 100 and 75 times faster when analysing 20 traits. The computational gains of HyPrColocterm·mth make it practical to perform multiple rounds of colocalizationterm·key analyses, each with different values of the prior p c and the threshold parameters P R * , P A * , to assess any sensitivity in the clusters of colocalized traits identified to changes in parameter specifications. An example of this, taken from data generated under scenario (iii), is presented in Fig. 7a . The resulting heatmap highlights the presence of four clusters of colocalized traits in the data and these clusters persist across most of the prior and threshold parameter settings. We include this sensitivity analysis in the HyPrColocterm·mth software and recommend its use. Fig. 7 HyPrColocterm·mth’s sensitivity analysis. Heatmap visualizing changes in the clusters of colocalized traits identified by HyPrColocterm·mth when using different choices of the colocalizationterm·key prior p c = {0.05, 0.02, 0.01, 0.005} and algorithm thresholds P R * = P A * = { 0.5 , 0.6 , 0.7 } . Cells appear darker when trait pairs cluster more often. Data were generated under scenario (iii) and when: a the single causal variantterm·key assumption is satisfied; or b the single causal variantterm·key assumption is violated.

We further tested the algorithm using a variety of thresholds { P R * , P A * } and two different prior frameworks (Supplementary Figs. S9 – S10 ). We also assessed results in the presence of correlated traits and overlapping samples ( Supplementary Information ). We analysed these data in three ways: (a) ignoring all correlation, i.e. wrongly assuming non-overlapping participants between pairs of studies and ignoring known trait correlation when setting the configuration prior probabilities; (b) adjusting for correlation between the summary data in the computation of the likelihood only; and (c) adjusting for correlation in the computation of the likelihood and accounting for known trait correlation when setting the configuration prior probabilities. Our findings suggest that analyses which account for correlation in the computation of the likelihood should also account for any known trait correlation in the configuration prior probabilities: the posterior probability of colocalizationterm·key between the truly colocalized traits in scenario (b), which ignored known correlation when setting the configuration prior, was significantly smaller than in scenario (c) - leading to a single large cluster of colocalized traits being split into smaller clusters (Supplementary Fig. S11 and Supplementary Table S2 ). Our results indicated that scenario (a), i.e. ignoring all correlation by treating studies as independent and traits as a-priori exchangeable, even when there is complete sample overlap (i.e. participants are the same in all studies), gives reasonable results and in our assessment was comparable to scenario (c) (Supplementary Fig. S10 and Tables S2 – S3 ). We discuss the theoretical reasons for this in the Supplementary Information . We additionally provide an example analysis protocol in our online vignette, which accompanies our software ( https://github.com/cnfoley/hyprcolocterm·mth ), offering further guidance on the choice of prior configuration probabilities and assessing any sensitivity of the clusters of colocalized traits identified to the choice of prior parameters.

Violations of the single causal variantterm·key assumption We assessed the performance of HyPrColocterm·mth when two or more traits have more than a single causal variantterm·key in the region. We simulated data for 10 traits and allowed up to 5 traits to have additional distinct causal variantsterm·key in the region, so that the sample contains a mixture of traits which either satisfy or violate the single causal variantterm·key assumption. The data are generated under three scenarios, as previously, but now each cluster of colocalized traits share a single causal variantterm·key and half of the traits in a cluster have secondary distinct causal variantsterm·key (‘Methods’). In terms of marginal genetic associations, the additional variants were randomly selected to explain either slightly less trait variance than the shared causal variantterm·key ( ≈ 0.75%) or the same amount of trait variance as the shared variant ( ≈ 1%).

The median accuracyterm·mth and TPRterm·mth of HyPrColocterm·mth reduced by as much as 38% - in scenario (i) - and had greater variation between the 1st and 9th deciles when the single causal variantterm·key assumption was violated (Supplementary Figs. S4a,b ); the reduction in performance was less pronounced when all studies had a large sample size. The FPRterm·mth remained modest however, i.e. the 1st decile and median FPRterm·mth were zero. A slight increase in the 9th decile of the FPRterm·mth was noted when causal variantsterm·key from distinct clusters were in strong LDterm·dom, i.e. r 2 > 0.95, removing these reduced the FPRterm·mth to zero (Supplementary Fig. S5c ). For larger samples sizes, the 1st, median and 9th deciles of the FPRterm·mth were approximately zero for each choice of prior (Supplementary Fig. S4c ). When considering only the clusters of traits identified as colocalizing with P R P A > 0.7, HyPrColocterm·mth again provided very reliable results across all three classification measures (Supplementary Fig. S6a–c ). Using the default settings { p = 10 −4 , p c = 0.02}, the algorithm generally performed well: in scenario (i) HyPrColocterm·mth regularly identified 8 of 10 traits as jointly colocalized; in scenario (ii) 5 out of 6 traits and; in scenario (iii) both clusters of colocalized traits, comprising 5 and 3 traits respectively (Supplementary Fig. S4f ) - highlighting HyPrColocterm·mth is conservative when additional causal variantsterm·key explain similar amounts of trait variation as the shared causal variantterm·key. We provide an illustration of HyPrColocterm·mth’s sensitivity analysis tool under scenario (iii) (Fig. 7b ) - correctly highlighting the presence of two clusters of colocalized traits. After applying more stringent prior and threshold values, one cluster reduced from 5 traits down to the 3 traits which have and share a single causal variantterm·key. This suggests strong evidence of 3 traits and weak evidence of 5 traits in the cluster. While the approach should be tailored to the problem at hand, if the analysis flags considerable sensitivity to the specification of the prior, we suggest: (a) reporting the clusters of colocalized traits identified as colocalizing with P R P A > 0.7 using the conservative prior setting p c = 0.02; and (b) where computationally practical, running pairwise analyses using a multi causal variantterm·key method, e.g. eCAVIAR 5 or ENLOC 6 , on the traits or clusters of traits which are reported in (a) but are not identified as colocalizing with P R P A > 0.7 using the more stringent prior p c = 0.01 - this may help clarify if traits are being removed from clusters owing to the presence of additional non-shared causal variantsterm·key, e.g. scenario (iii) (Fig. 7b ), and should therefore be reported. We provide further guidance on the reliability of the BB algorithm when secondary causal variantsterm·key are added to all traits in the region and when varying LDterm·dom between causal variantsterm·key ( Supplementary Information ; Supplementary Table S5 ).

We also compared results with those obtained using pairwise COLOC and eCAVIAR 5 (with a colocalizationterm·key posterior probability, CLPP, cut-off of 1% and default prior choices), another software package for colocalizationterm·key which allows each trait to have multiple causal variantsterm·key but is limited to the analysis of pair of traits only. We note that the SNP level CLPP measure of eCAVIAR is computed in the presence of multiple causal variantsterm·key and is distinct from the SNP level probabilities, computed under a single causal variantterm·key assumption, which comprise the posterior probability measure returned by HyPrColocterm·mth and COLOC - making comparisons between the methods challenging. We compare the methods as they are used in practice, summarizing HyPrColocterm·mth and COLOC using the posterior probability of the hypothesis that a cluster or a pair of traits colocalize 2 , 8 , 14 and summarizing eCAVIAR using the SNP-level CLPP. Our choice of CLPP cut-off of 1% was shown to have a low FPRterm·mth across a range of scenarios previously 5 . In our analyses we found that pairwise eCAVIAR had increased accuracyterm·mth relative to HyPrColocterm·mth and pairwise COLOC, e.g. in scenario (i) median accuracyterm·mth improved by as much as 0.15 (when sample sizes varied) and 0.2 (when sample sizes were large) (Supplementary Fig. S4a and Table S8 ). Broadly, this was a result of the single causal variantterm·key methods having a lower TPRterm·mth (Supplementary Fig. S4a, b ). However, by borrowing information between multiple traits HyPrColocterm·mth outperformed eCAVIAR when fine-mappingterm·key the shared causal variantterm·key (Supplementary Fig. S4d )—despite not incorporating LDterm·dom information. After thresholding the posterior to P R P A > 0.7, HyPrColocterm·mth again outperformed pairwise COLOC (Supplementary Fig. S6a–c ).

Despite violations of the single causal variantterm·key assumption, our analyses demonstrate that HyPrColocterm·mth can continue to identify clusters of colocalized traits, returning conservative results otherwise, with major computational advantages over competing software: in the analysis of 10 traits and in a region containing around 1000 SNPs, the single joint colocalizationterm·key analysis of HyPrColocterm·mth was computed approximately 100,000 times faster than the 45 pairwise analyses of eCAVIAR. The HyPrColocterm·mth algorithm can additionally be used to rapidly identify genomic regions and clusters of traits to better prioritize the use of more computationally expensive multi-causal variant colocalizationterm·key software for pairs of traits ( Supplementary Information ).

Map of genetic risk shared across CHDterm·dom and related traits We used HyPrColocterm·mth to investigate genetic associations shared across CHDterm·dom 17 and 14 related traits: 12 CHDterm·dom risk factors 18 – 22 , a comorbidity 23 and a social factor 24 (Supplementary Table S1 for details). We performed colocalizationterm·key analyses in pre-defined disjoint LDterm·dom blocks spanning the entire genome 25 . To highlight that multi-trait colocalizationterm·key analyses can aid discovery of new disease-associated loci, we used the CARDIoGRAMplusC4D 2015 data for CHDterm·dom 17 , which brought the total number of CHD-associated regions to 58, and contrasted our findings with the current total of ~160 CHD-associated regions 26 . For each region in which CHDterm·dom and at least one related trait colocalized, we integrated whole blood gene expression 27 quantitative trait loci (eQTL) and protein expression 28 quantitative trait loci (pQTL) information into our analyses to prioritise candidate causal genes (‘Methods’).

Multi-trait colocalizationterm·key Our genome-wide analysis identified 43 regions in which CHDterm·dom colocalized with one or more related traits (Fig. 8 and Tables 1 – 3 ). Twenty-three of the 43 colocalizationsterm·key involved blood pressure, consistent with blood pressure being an important risk factor for CHDterm·dom 29 . Other traits colocalizing with CHDterm·dom across multiple genomic regions were cholesterol measures (16 regions); adiposity measures (9 regions); type 2 diabetes (T2D; 4 regions) and; rheumatoid arthritis (2 regions). Moreover, by colocalizing CHDterm·dom and related traits, our analyses suggest these traits share some biological pathways. Fig. 8 Genome-wide multi-trait colocalizationterm·key analysis of CHDterm·dom and fourteen related traits. a Summary of the number of regions across the genome in which CHDterm·dom colocalizes with at least one related trait. Results are aggregated by trait family, e.g. major lipids, and by each individual trait (see Supplementary Table S1 for a list of trait abbreviations). b Stacked association plots of CHDterm·dom with high-density lipoprotein (HDL), low-density lipoprotein (LDL), systolic blood pressure (SBP), diastolic blood pressure (DBP) and rheumatoid arthritis (RA). HyPrColocterm·mth implicated both the SH2B3-ATXN2 locus and risk variant rs713782, both of which have been previously reported as associated with CHDterm·dom risk 26 . However, HyPrColocterm·mth extended this result by identifying that the risk loci and variant are shared with 5 conventional CHDterm·dom risk factors 11 . SNPs in stronger LDterm·dom with the putative causal SNP rs713782 appear darker in the plot. c HyPrColocterm·mth identified rs713782 as a candidate causal variantterm·key explaining the shared association signal between CHDterm·dom and the 5 related traits. The posterior probability of colocalizationterm·key between the traits was 0.909 and rs713782 explained over 76% of this, i.e. the posterior probability of rs713782 being the shared causal variantterm·key is 0.909 × 0.76 = 0.69. The next candidate variant explained <20%. Table 1 CHDterm·dom loci that were known at the time of the CARDIoGRAMplusC4D data release (2015). Chr Locus Traits Colocalized SNP (consequence) Gene Known CHDterm·dom locus (SNP) PPFC (PPE) Expressed gene (eQTL) Protein (pQTL) 2 ABCG8, ABCG5 CHDterm·dom, LDL rs4299376 (Intron) ABCG8 Yes 32 (Yes 32 ) 0.918 (0.949) – – 4 GUCY1A1 CHDterm·dom, DBP rs72689147 (Intron) GUCY1A1 Yes 16 (Yes 17 ) 0.931 (0.241) GUCY1A1 (rs12643599) – 6 PHACTR1, EDN1 CHDterm·dom, SBP rs9349379 (Intron) PHACTR1 Yes 33 , 35 (Yes 33 ) 0.999 (1) – – 6 LPA CHDterm·dom, LDL rs10455872 (Intron) LPA Yes 32 , 35 (Yes 32 , 35 ) 0.998 (0.538) – – 7 HDAC9 CHDterm·dom, SBP rs2107595 (Intergenic) HDAC9 Yes 16 (Yes 17 ) 0.996 (0.729) – – 7 ZC3HC1, KLHDC10 CHDterm·dom, DBP rs11556924 (Missense) ZC3HC1 Yes 16 , 32 , 35 (Yes 16 , 32 , 35 ) 0.999 (0.994) – – 8 TRIB1 CHDterm·dom, HDL, LDL, TG, eGFR rs2954029 (Intron) RP11-136O12.2 Yes 16 (Yes 16 ) 0.925 (0.872) – – 9 ANRIL, CDKN2B-AS1 CHDterm·dom, DBP rs2891168 (Intron) CDKN2B-AS1 Yes 17 (Yes 17 ) 0.870 (0.755) – – 9 ABO CHDterm·dom, LDL, DBP, T2D rs507666 (Intron) ABO Yes 16 , 35 (Yes 17 ) 0.984 (0.582) – – 10 KIAA1462 CHDterm·dom, DBP rs1887318 (Intron) KIAA1462 Yes 16 , 33 (Yes 17 ) 0.937 (0.433) – – 11 APOA1-C3-A4-A5 CHDterm·dom, HDL, LDL, TG rs964184 (3 prime UTR) ZPR1, BUD13 Yes 35 (Yes 35 ) 0.957 (1) – Apolipo-protein A-V (rs964184) 12 ATP2B1 CHDterm·dom, SBP rs2681492 (Intron) ATP2B1 Yes 17 (Yes 17 ) 0.980 (0.303) – – 12 SH2B3 CHDterm·dom, HDL, LDL, SBP, DBP, RA rs7137828 (Intron) ATXN2 Yes 35 (Yes 17 ) 0.909 (0.768) TRAFD1 (rs7137828) – 15 FES, FURIN CHDterm·dom, SBP, DBP rs35346340 (Splice region) FES Yes 16 (Yes 17 ) 0.959 (0.579) FES (rs8027450) – 18 MC4R, PMAIP1 CHDterm·dom, HDL, TG, BMI, WC rs12967135 (Intergenic) – Yes 17 (Yes 17 ) 0.859 (0.434) – – 19 LDLR, SMARCA4 CHDterm·dom, LDL rs112374545 (Intergenic) LDLR Yes 16 , 35 (Yes 17 ) 0.937 (0.556) – – 19 APOC1, APOE, PVRL2, COTL1 CHDterm·dom, HDL, WC rs4420638 (Downstream) APOC1 Yes 17 (Yes 17 ) 0.959 (0.999) – Apolipo-protein E (rs4420638) 21 KCNE2 CHDterm·dom, DBP rs28451064 (Intron) AP000318.2 Yes 17 (Yes 17 ) 0.998 (0.974) – – HyPrColocterm·mth identified eighteen known CHDterm·dom genetic risk loci (i.e. CHDterm·dom loci reported before or at the time of the CARDIoGRAMplusC4D data release in 2015) with colocalized associations across CHDterm·dom and one or more of 14 related traits. Chr: denotes chromosome; Locus: candidate causal gene(s) as listed by Erdmann et al. 40 ; Traits: traits with colocalized association; Colocalized SNP(consequence): SNP marking association shared across the traits and its annotation in VEP 50 from PhenoScanner 51 ; Gene: nearest gene to colocalized SNP; Known CHDterm·dom locus: locus known at time of 2015 CHDterm·dom data release 17 (i.e. published in ref. 17 or earlier) or subsequently identified 26 ; PPFC: posterior probability of colocalizationterm·key; PPE: proportion of PPFC explained by the listed SNP; eQTL: gene expression 27 ; pQTL: protein expression 28 . See Supplementary Table S1 for a list of the trait abbreviations. Full results from these analyses are available at https://jrs95.shinyapps.io/hyprcolocterm·mth_chdterm·dom . Table 2 CHDterm·dom loci reported after the time of the CARDIoGRAMplusC4D data release (2015). Chr Locus Traits Colocalized SNP (consequence) Gene Known CHDterm·dom locus (SNP) PPFC (PPE) Expressed gene (eQTL) Protein (pQTL) 1 PRDM16 CHDterm·dom, SBP, DBP rs2493288 (Intron) PRDM16 Yes 26 (Yes 26 ) 0.8009 (0.3471) – – 1 FHL3 CHDterm·dom, SBP rs34655914 (Missense) INPP5B Yes 26 (Yes 26 ) 0.9468 (0.0832) SF3A3 (rs28428561) ; UTP11L (rs4360494) ; RNU6-510P (rs61776719) – 1 SORT1 CHDterm·dom, HDL rs12740374 (3 prime UTR) CELSR2 Yes 26 (Yes 26 ) 0.9898 (0.9997) – – 1 LMOD1 CHDterm·dom, BMI, WC rs2678204 (Intron) IPO9 Yes 30 (Yes 30 ) 0.8273 (0.1627) IPO9 (rs2494115) – 2 FIGN CHDterm·dom, SBP rs268263 (Intron) AC092684.1 Yes 26 (Yes 26 ) 0.789 (0.995) – – 2 IRS1 CHDterm·dom, HDL, TG rs62188784 (Intergenic) AC068138.1 Yes 26 (Yes 26 ) 0.8234 (0.4852) – – 3 RHOA CHDterm·dom, BMI, EDU rs73078367 (Downstream) NCKIPSD Yes 26 (Yes 26 ) 0.9541 (0.5656) – – 3 RHOA CHDterm·dom, SBP rs7623687 (Intron) RHOA Yes 34 (Yes 34 ) 0.9713 (0.2455) – – 4 FGF5, PRDM8 CHDterm·dom, SBP, DBP rs13125101 (Intergenic) FGF5 Yes 26 (Yes 26 ) 0.9827 (0.4148) – – 5 MAP3K1 CHDterm·dom, HDL, TG, WC, SBP, T2D rs9686661 (Intron) C5orf67 Yes 26 (Yes 26 ) 0.7755 (0.7115) – – 6 VEGFA CHDterm·dom, HDL, TG, BMI, WC rs998584 (Downstream) VEGFA Yes 26 (Yes 26 ) 0.8376 (0.9746) – – 10 TSPAN14, FAM213A CHDterm·dom, RA rs2343306 (Intron) TSPAN14 Yes 26 (No) 0.9064 (0.7279) – – 11 ARNTL CHDterm·dom, DBP rs10832013 (Upstream) ARNTL Yes 26 (Yes 26 ) 0.9403 (0.0823) – – 11 SIPA1 CHDterm·dom, HDL, TG rs12801636 (Intron) PCNX3 Yes 30 (Yes 30 ) 0.8369 (0.8945) – – 12 HNF1A CHDterm·dom, LDL rs1169288 (Missense) HNF1A Yes 30 (Yes 30 ) 0.9645 (0.5762) – – 13 N4BP2L2, PDS5B CHDterm·dom, BMI rs35193668 (Intron) N4BP2L2 Yes 26 (Yes 26 ) 0.6785 (0.0911) N4BP2L2 (rs9337) – 16 CDH13 CHDterm·dom, DBP rs7500448 (Intron) CDH13 Yes 26 (Yes 26 ) 0.9947 (1) – – 16 CTRB2, BCAR1 CHDterm·dom, T2D rs55993634 (Downstream) CTRB2 Yes 34 (Yes 26 ) 0.8296 (0.3868) BCAR1 (rs28595463) – 17 IGF2BP1 CHDterm·dom, BMI, T2D rs11079849 (Intron) IGF2BP1 Yes 26 (Yes 26 ) 0.8389 (0.831) – – 17 PECAM1, DDX5, TEX2 CHDterm·dom, SBP, DBP rs1867624 (Upstream) RPL31P57 Yes 30 (Yes 30 ) 0.7963 (0.4276) - – HyPrColocterm·mth identified twenty CHDterm·dom genetic risk loci—reported after the time of the CARDIoGRAMplusC4D data release in 2015—with colocalized associations across CHDterm·dom and one or more of 14 related traits. See Table 1 for a full description of the table items. Table 3 New CHDterm·dom loci sharing colocalized associations with related traits. Chr Locus Traits Colocalized SNP (consequence) Gene Known CHDterm·dom locus (SNP) PPFC (PPE) Expressed gene (eQTL) Protein (pQTL) 6 FHL5 CHDterm·dom, SBP rs9486719 (Intron) FHL5 – 0.844 (0.1542) – – 10 CYP26A1 CHDterm·dom, TG rs2068888 (Downstream) CYP26A1 – 0.8454 (0.7669) – – 16 ANKRD11 CHDterm·dom, WC rs11643561 (Intron) ANKRD11 – 0.7827 (0.0795) – – 19 RSPH6A CHDterm·dom, SBP rs8108474 (Intron) RSPH6A – 0.7802 (0.1435) – – 20 PREX1 CHDterm·dom, SBP, DBP rs79044887 (Intron) PREX1 – 0.7237 (0.132) – – HyPrColocterm·mth identified five regions—not yet reported as CHDterm·dom genetic risk loci—with colocalized associations across CHDterm·dom and one or more related trait. See Table 1 for a full description of table items.

In thirty-eight of the 43 (88%) colocalized regions, the candidate causal SNP proposed by HyPrColocterm·mth and/or its nearest gene, have been previously identified 16 , 17 , 26 , 30 – 35 . Importantly, 20 of these were reported after the CARDIoGRAMplusC4D study 17 . For example, FGF5 was sub-genome-wide significant ( P > 5 × 10 −8 ) with CHDterm·dom in the 2015 data, but through colocalizationterm·key with blood pressure, we highlight it as a CHDterm·dom locus and it is genome-wide significant in the most recent CHDterm·dom GWASterm·mth 26 . The remaining 18 regions were reported previously, but one, APOA1-C3-A4-A5 , was sub-genome-wide significant in the CARDIoGRAMplusC4D study 17 despite having been reported previously 35 . However, we used HyPrColocterm·mth to show that the association of major lipids colocalize with a CHDterm·dom signal, highlighting this as a CHDterm·dom locus in these data (Table 1 and Supplementary Fig. S13 ). The locus has subsequently been replicated 26 , 31 and we show below that the signal also colocalizes with circulating apolipoprotein A-V protein levels (Table 1 ). This demonstrates that joint colocalizationterm·key analyses of diseases and related traits can improve power to detect new associations (an approach which is advocated outside of colocalizationterm·key studies 36 ). Our results also illustrate that multi-trait colocalizationterm·key analyses can provide further insights into well-known risk-loci of complex disease. For example, at the well-studied SH2B3 - ATXN2 region 26 , 35 , HyPrColocterm·mth detected two cholesterol measures (LDL, HDL), two blood pressure measures (SBP, DBP) and rheumatoid arthritis (RA) colocalizing with CHDterm·dom at the previously reported CHD-associated SNP 26 rs7137828 (PPFC = 0.909 of which 76.8% is explained by the variant rs7137828; Fig. 8 ). In addition, we implicated a candidate SNP and locus in a further 5 CHDterm·dom regions not previously associated with CHDterm·dom risk (Table 3 ). In one of the 5 regions, CYP26A1 , CHDterm·dom colocalized with tri-glycerides (TG) and HyPrColocterm·mth identified a single variant that explained over 75% of the posterior probability of colocalizationterm·key, supporting this SNP as a candidate shared CHDterm·dom/TG variant.

For each of the 43 regions that shared genetic associations across CHDterm·dom and related traits, we further integrated whole blood gene 27 and protein 28 expression into the colocalizationterm·key analyses. We tested cis eQTL for 1828 genes and cis pQTL from the 854 published proteins across the 43 loci for colocalizationterm·key with CHDterm·dom and the related traits. Of the 43 listed variants (Tables 1 – 3 ), 27 were associated with expression of at least one gene ( P < 5 × 10 −8 ) and a total of 125 such genes were identified. HyPrColocterm·mth refined this, identifying six regions colocalizing with eQTL for one expressed gene and one region, the FHL3 locus, colocalizing with expression of three genes ( SF3A3 , UTP11L , RNU6-510P ) (Table 2 ). The GUCY1A3 locus has previously been associated with BP 37 and with CHDterm·dom 16 . Here we show that these associations are likely to be due to the same variant, rs72689147 (PPFC = 0.93), with the G allele increasing DBP and risk of CHDterm·dom. We furthermore show that the association colocalizes with expression of GUCY1A1 in whole blood, with the G allele reducing GUCY1A1 expression (PPFC = 0.77; Table 1 ). The GUCY1A1 gene is ubiquitously expressed in heart tissues, including in the coronary and aortic arteries 38 . In the mouse, higher expression of GUCY1A1 has been correlated with less atherosclerosis in the aorta 39 . GUCY1A1 is a likely candidate gene in this locus 40 , illustrating the utility of HyPrColocterm·mth to help prioritise candidate causal genes. The CTRB2 - BCAR1 locus was not known at the time of the release of the 2015 CARDIoGRAMplusC4D data, however we find the association at this locus is shared with T2D (PPFC = 0.83) and that BCAR1 expression colocalized with the CHDterm·dom association (PPFC = 0.86). Other studies have implicated the locus in CHDterm·dom 34 and suggested BCAR1 as the causal gene in carotid intimal thickening 41 , 42 . We note that two CHDterm·dom loci also colocalize with circulating plasma proteins, APOA1-C3-A4-A5 , with apolipoprotein A-V and the APOE locus with apolipoprotein E (Table 1 ).

Of the 38 known CHDterm·dom loci that colocalized with a related trait, 8 are reported to have a single causal variantterm·key 26 , of these we identified the same CHD-associated variant (or one in LDterm·dom with either r 2 > 0.8 or |D′|> 0.8) 15 at seven loci ( SORT1 , PHACTR1 , ZC3HC1 , CDKN2B-AS1 , KCNE2 , CDH13, APOE ). Despite the possible presence of multiple causal associations at other loci, HyPrColocterm·mth was still able to pick out single shared associations across traits: a result supported by our simulation study when additional distinct causal variantsterm·key explain less or similar trait variation than that explained by a shared causal variantterm·key between colocalized traits ( Supplementary Information ). In our analyses, we set p c = 0.02, i.e. γ = 0.98, and report only the clusters of traits whose posterior probability of colocalizationterm·key was greater than 0.7. We assessed sensitivity to the choice of colocalizationterm·key prior, repeating analyses with p c = 0.01, and found no appreciable difference in the clusters identified (results not reported).

Discussion We have developed and applied a deterministic Bayesianterm·key colocalizationterm·key algorithm, HyPrColocterm·mth, for multi-trait statistical colocalizationterm·mth analyses. HyPrColocterm·mth is based on the same underlying statistical model as COLOC 2 , but enables colocalizationterm·key analyses to be performed across massive numbers of traits, owing to the insight that the posterior probability of colocalizationterm·key at a single causal variantterm·key can be accurately approximated by enumerating only a small number of putative causal configurations. HyPrColocterm·mth avoids repeated rounds of pairwise colocalizationterm·key analyses which can inflate the false negative rate and have reduced performance in identifying a shared causal variantterm·key. The HyPrColocterm·mth algorithm was validated using simulations and used to assess genetic risk shared across CHDterm·dom and related traits. Using CHDterm·dom data from 2015 17 , in which 46 regions were genome-wide significant ( P < 5 × 10 −8 ), our multi-trait colocalizationterm·key analysis identified 43 regions in which CHDterm·dom colocalized with ≥1 related trait. With this approach, we were able to identify CHDterm·dom loci that were not known at the time of the data release (2015), demonstrating the benefit of synthesising data on related traits to uncover potential new disease-associated loci 8 , 36 . A further five regions, we postulate, may be identified as CHDterm·dom loci in the future. Others have considered pleiotropic effects of CHDterm·dom loci previously 43 , but our formal colocalizationterm·key analyses are more robust, e.g. in the ABO region we show colocalizationterm·key of T2D and DBP in addition to the previously reported pleiotropic effect with LDL. We integrated eQTL and pQTL data to prioritise candidate genes at some loci, e.g. GUCY1A1 , BCAR1 and APOE .

The HyPrColocterm·mth algorithm identifies regions of the genome where there is evidence of a shared causal variantterm·key (by dissecting the genome into distinct regions) and also allows for a targeted analysis of a specific genomic locus of primary interest, e.g. when aiming to identify the perturbation of a biological pathway through the influence of a particular gene. Moreover, these region-specific analyses can highlight candidate causal genes, which will help improve biological understanding and may indicate potential drug targets to inform medicines development 44 .

We have described HyPrColocterm·mth under the assumption of at most one causal variantterm·key per trait. Future work is required to extend this methodology and algorithm to multiple-causal variants. We note that the reliability of results under the single causal variantterm·key assumption only break down when secondary causal variantsterm·key explain as much trait variation as the shared variant ( Supplementary Information ). An example of which is the expression of SH2B3 , where multiple causal variantsterm·key for the expression of this gene masks colocalizationterm·key with the CHDterm·dom signal, we discuss an approach to building colocalizationterm·key analyses which might help support the single causal variantterm·key assumption ( Supplementary Information ). We note that misspecification of LDterm·dom between causal variantsterm·key has a major impact on correct detection of multiple causal variantsterm·key in a region 45 , making a single causal variantterm·key assessment the most reliable when accurate study-level LDterm·dom information is not available. To overcome challenges when specifying the prior probability of a causal configuration, we have suggested two different parsimonious configuration priors (‘Methods’). The computational advantages of HyPrColocterm·mth make it practical to assess sensitivity of results to the specification of prior and threshold parameters as part of regular use. The HyPrColocterm·mth software includes a tool to do this, visualizing any changes in the clusters of colocalized traits identified as parameters are varied. Nevertheless, other priors may be more appropriate for particular applications.

In summary, we have developed a computationally efficient method that can perform multi-trait colocalizationterm·key on a large scale. As the size and scale of available data on genetic associations with traits increase, computationally scalable methods such as HyPrColocterm·mth will be increasingly valuable in prioritizing causal genes and revealing causal pathways.

Methods SNP association models Let Y i denote one of i = 1,2,…, m , traits assessed in a maximum of m studies, i.e. two or more traits can be measured in the same study, and G ij denote the genotype of the j th genetic variant. It is assumed that the outcome model for Y i is given by 4 E Y i G i j = h i − 1 α i j + β i j G i j , where α ij is the intercept term and h i is a function linking the i th outcome to the genotype G ij , for all j = 1,2,…, Q genetic variants in the genomic region. The function h i is typically taken as the identity function for continuous traits and the logit function for binary traits. The aim of colocalizationterm·key analyses is to identify genomic loci where there exists an G ij that is causally associated with at least two of the m traits. For each of the m traits and Q genetic variants, we assume that GWAS summary statisticsterm·mth β ^ i j and var β ^ i j are available. We use these data to perform colocalizationterm·key analyses in genomic loci.

Colocalizationterm·key posterior probability Using binary vectors to indicate whether a variant putatively causally influences a trait, we can define causal configurations ( S ) that can be grouped into sets ( S H ) which belong to a single data generating hypothesis ( H ). We use the notation H i , j , … to denote a set of hypotheses in which a collection of i traits share a causal variantterm·key, a separate collection of j traits share a distinct causal variantterm·key, and so on (Fig. 1 ). For, example, H 2 , 1 denotes the set of hypotheses in which each hypothesis specifies uniquely 2 traits that share a causal variantterm·key, a single trait has a distinct causal variantterm·key and all remaining m −3 traits do not have a causal variantterm·key in the region. Assuming at most one causal variantterm·key for each trait these data generating hypotheses can be combined to generate a hypothesis space (Ω). The posterior probability of hypothesis H , given the combined data D from all m studies, can therefore be computed using ( Supplementary Information ), 5 P H ∣ D = ∑ S ∈ S H B F ( S ) p ( S ) p ( S 0 ) ∑ H i ∈ Ω ∑ S ∈ S H i B F ( S ) p ( S ) p ( S 0 ) , where p ( S )/ p ( S 0 ) is the prior-odds of configuration S ∈ S H compared with the null-configuration S 0 , i.e. no genetic association with any trait. See ref. 2 for a derivation with m = 2 traits. BF ( S ) is a Bayes factor which is the likelihood of the data being generated under S ∈ S H relative to the likelihood of the data being generated S 0 .

We describe the space of multi-trait colocalizationterm·key models using a set of mutually exclusive hypotheses and causal configurations as this approach extends the methodology and language used previously 2 , 8 . We note, however, that each causal configuration is equivalent to a model which, for each trait, details the location of the causal variantterm·key in the region. Hence, the problem of identifying a hypothesis and causal configuration with the greatest support given the data D , is equivalent to identifying the joint trait-variant model with greatest support 2 , 13 .

Computing Bayes factors: independent studies If the trait associations are calculated using independent studies (i.e. no overlapping samples in the GWASterm·mth datasets), the Bayes factors can be computed using Wakefield’s Approximate Bayes Factors 13 ( ABF ) for each trait i and genetic variant j , i.e. 6 A B F i j = v i j 2 v i j 2 + w i j 2 exp z i j 2 2 × w i j 2 v i j 2 + w i j 2 , where z ij , v ij and w ij are the Z-statistic, standard error and the prior standard deviation for β ^ i j , respectively. Following 2 , for continuous variables w ij is set to 0.15 while for binary traits it is set to 0.2. As an example, the ABF for the hypothesis that all m traits colocalize at genetic variant j ( S j ∈ S m ) is given by, 7 A B F S j = ∏ i m A B F i j .

Calculating Bayes factors: non-independent studies If the trait associations are not calculated using independent studies i.e. there are overlapping samples, the Bayes factor for each causal configuration can be computed using a Joint ABF ( JABF ) ( Supplementary Information ). The JABF for causal configuration S is given by, 8 J A B F S = Σ β ^ Σ β ^ + Σ β exp 1 2 β ^ T Σ β ^ + Σ β − 1 Σ β Σ β ^ − 1 β ^ , where β ^ is the vector of regression coefficients for all m traits, Σ β ^ is an m × m variance-covariance matrix of the regression coefficients (i.e. V ρ ^ V , where V 2 is a diagonal matrix of variances for the regression coefficients, e.g. with i th diagonal element v i ⋅ 2 , and ρ ^ is the observed correlation matrix for the regression coefficients) and Σ β is the prior variance-covariance matrix (i.e. W ρ W , where W 2 is a diagonal matrix of prior variance, e.g. with i th diagonal element w i ⋅ 2 , and ρ is the prior correlation matrix between traits). The correlation matrix ( ρ ^ ) is computed using the tetrachoric correlation method 46 and we discuss our approach to setting ρ in the Supplementary Information .

Configuration prior probabilities We consider two different strategies for determining the priors for different hypotheses: variant-level priors and uniform priors.

Variant-level prior probabilities The prior probability space for a single genetic variant can be fully partitioned into the prior probability that the genetic variant is not associated with any of the m traits, p 0 , the prior probability that the genetic variant is associated with only the first trait, p 1 ,…, the prior probability that the SNP is associated with a subset of k traits { j 1 , j 2 ,…, j k }, p j 1 j 2 … j k , …, the prior probability that the genetic variant is associated with all traits, p 12… m . Hence, 9 p 0 + ∑ k = 1 m ∑ j 1 = 1 m ∑ j 2 > j 1 … ∑ j k > j k − 1 p j 1 j 2 … j k = 1 .

The space therefore requires the specification of 2 m prior parameters which, even for modest values of m , is computationally impractical. Following 2 , 8 we set that the prior probability to not vary by genetic variant, nor by the specific collection of colocalized traits of a given size, but by the number of colocalized traits, i.e. a SNP associated with a total of k traits has a prior probability that depends on the number k but not the specific collection of traits. To allow for the assessment of large numbers of traits we propose variant-level priors where the prior probability that a genetic variant is associated with k traits is given by, 10 p 12 … k = p ∏ i = 2 k 1 − γ i − 1 , k = 2 , … m , where p is the probability of the genetic variant being associated with one trait and γ is a parameter which controls the probability that a genetic variant is associated with an additional trait. Notably, 1 − γ is the probability of a variant being causal for a second trait given it is causal for one trait, i.e. it is the conditional colocalizationterm·key prior p c , p c = 1 − γ , 1–γ 2 is the probability it is causal for a third trait given it is causal for two traits, and so on.

It follows that, 11 p S p S 0 = p 12 … k p 0 = p p 0 ∏ i = 2 k 1 − γ i − 1 , k = 2 , … , m , for configurations S ∈ S H k , where k traits share a causal variantterm·key and the remaining m − k traits do not have a casual variant, and 12 p S p S 0 = p 12 … m − 1 p 1 p 0 2 = p p 0 2 ∏ i = 2 m − 1 1 − γ i − 1 , for configurations S ∈ S H ( m − 1 , 1 ) , where m − 1 traits share a causal variantterm·key and the remaining trait has a distinct causal variantterm·key. This prior set-up allows evidence to grow in favour of k traits colocalizing conditional on evidence supporting k − 1 traits colocalizing ( Supplementary Information ). For example, if the first k traits are believed to share a causal variantterm·key a-priori, then the prior probability that the ( k + 1) th is also colocalized, conditional on the other k traits, increases as the number of colocalized traits k grows. The marginal prior probability of k traits colocalizing is always very small, however, which controls the false positive rateterm·mth (Fig. 6 and Supplementary Figs. S2 – 6 ; Supplementary Tables S2 – 3 ). Conditional growth limits the loss of power when assessing colocalizationterm·key across a large number of traits. A loss in power necessarily occurs when analysing large numbers of colocalized traits, due to the rapid growth in the number of hypotheses in which a subset of traits can colocalize relative to all traits colocalizing. Evidence supporting these ‘subset’ hypotheses will eventually overwhelm evidence in favour of the maximum number of truly colocalized traits for a fixed sample size (top row Fig. 4 ). Based on our simulation results (Fig. 6 and Supplementary Figs. S2 – 6 ) and previous investigations 14 , we recommend users set p c = 0.02, i.e. γ = 0.98, and report results from the clusters of colocalized traits which satisfy P R P A > 0.7. Setting p c = 0.02 increased the number of datasets in which clusters satisfying P R P A > 0.7 were identified (c.f. simulation study) while maintaining a low FPRterm·mth throughout. Using the same posterior threshold of 0.7 and setting p c = 0.05 returned reasonable results. However, we do not recommend users set p c = 0.05 due to the slight increase in the 9th decile of the FPRterm·mth in scenario (iii) (Fig. 6c ). If two or more traits in a cluster are known to be related, this information would ideally be included in analyses and we outline an extension to our prior setup which allows for non-exchangeability of traits to be included ( Supplementary Information ).

Conditionally uniform prior probabilities An alternative prior strategy is to assume uniform priors for each configuration within a hypothesis 47 . This strategy benefits from: (i) not setting variant-level information and (ii) implicitly accounting for large differences in the causal configuration space between hypotheses, which limits the loss in power of the PPFC for very large m . These priors take the form, 13 P S ∣ H P S 0 ∣ H 0 = 1 / ∣ S H ∣ 1 / ∣ S 0 ∣ = 1 / ∣ S H ∣ , where S H k = Q and 14 S H ( m − 1 , 1 ) = Q Q − 1 : m = 2 , m Q Q − 1 : m > 2 .

Through simulations, we identified the conditionally uniform prior as less conservative than variant-level priors, having an increased false detection rate of colocalizationterm·key. ( Supplementary Information ; Supplementary Figs. S10 , 11 ). This could lead to an increased false positive detection rate in practice.

HyPrColocterm·mth posterior approximation To compute the posterior probability of full colocalizationterm·mth across a large number of traits we propose the HyPrColocterm·mth posterior approximation. Let P ( H m | D ), P scv , P ( m −1,1) and P all denote: (i) the posterior probability of full colocalizationterm·mth; (ii) the sum of the posterior probabilities in which no traits have a causal variantterm·key, a subset of m − 1 traits share a causal variantterm·key (the remaining trait does not have a causal variantterm·key) and all m traits colocalize ( P scv ); (iii) the sum of posterior probabilities in which a subset of m − 1 traits share a causal variantterm·key and the remaining trait has a distinct causal variantterm·key ( P ( m −1,1) ) and; (iv) the sum of all posterior probabilities of at most one causal variantterm·key per trait ( P all ). That is, 15 P s c v = P H 0 ∣ D + P ( H m − 1 ∣ D ) + P H m ∣ D and P ( m − 1 , 1 ) = P ( H ( m − 1 , 1 ) ∣ D ) .

The HyPrColocterm·mth posterior is computed in two steps. Step 1 computes the regional association probability P R , defined as: 16 P R = P ( H m ∣ D ) P s c v ≥ P ( H m ∣ D ) .

Step 2 computes the alignment probabilityterm·mth P A , defined as: 17 P A = P ( H m ∣ D ) P ( H m ∣ D ) + P ( m − 1 , 1 ) ≥ P H m ∣ D .

Note that P R is computed using ( m + 1) Q causal configurations and P A is computed using an additional mQ ( Q − 1) causal configurations. Hence, computation of P R and P A has O m Q 2 computational cost. We let P a l l c = P a l l − P s c v − P m − 1 , 1 , then it follows that the posterior probability of all traits sharing a single causal variantterm·key is given by 18 P H m ∣ D = P H m ∣ D P a l l = P ( H m ∣ D ) P s c v P s c v P a l l = P ( H m ∣ D ) P s c v × P s c v P ( H m ∣ D ) P ( H m ∣ D ) P s c v P ( H m ∣ D ) P ( H m ∣ D ) + P ( m − 1 , 1 ) − P s c v P ( H m ∣ D ) 1 − P ( H m ∣ D ) P s c v P ( m − 1 , 1 ) − P ( H m ∣ D ) P s c v P a l l c = P R P A 1 − 1 − P R 1 − P A − P R 1 − P A P a l l c P m − 1 , 1 = P R P A + O δ A 2 + δ R δ A , δ R , δ A → 0 , where δ R = 1− P R , δ A = 1− P A and P a l l c P ( m − 1 , 1 ) = O δ R + δ A , ( Supplementary Information ). By definition, P ( H m ∣ D ) → 1 ⇔ P R → 1 and P A → 1 . Hence together the regional and alignment probabilities when multiplied form a statistic that is sufficient to accurately assess evidence of the full colocalizationterm·key hypothesis. The objects P R and P A can be defined for various collections of hypotheses that partition P all . However, the major insight is that the hypotheses contained in P R and P A are computed with minimal computation burden, i.e. computed using ≤ mQ 2 causal configurations, amongst all alternatives, making the HyPrColocterm·mth approximation tractable for very large numbers of traits m .

Our software allows for the assessment of the HyPrColocterm·mth approximation by increasing the number of hypotheses used to approximate P R , e.g. we can compute 19 P R ′ = P ( H m ∣ D ) P H 0 ∣ D + P H m − 2 ∣ D + P H m − 1 ∣ D + P ( H m ∣ D ) , which is computed from O m 2 Q causal configurations and assess the relative difference between P R and P R ′ . We show that P R ′ = P R 1 + δ R ( Supplementary Information ) and through simulations that there very close correspondence between P R ′ and P R (Supplementary Table S4 ).

Branch and Bound divisive clustering algorithm To identify complex patterns of colocalizationterm·key amongst all traits, we propose a branch and bound (BB) divisive clustering algorithm that utilizes the HyPrColocterm·mth approximation to identify a cluster of traits with the greatest evidence of colocalizationterm·key at each iteration (Supplementary Fig. S1a and Supplementary Information ). Starting with all of the traits in a single cluster, the algorithm explores evidence supporting any of 2 m branches - a branch represents a hypothesis whereby m − 1 traits share a causal variantterm·key and either the remaining trait does not have a causal variantterm·key or has a causal variantterm·key elsewhere in the region - against the full colocalizationterm·key hypothesis. These branches represent the hypotheses used in the computation of the regional and alignment probabilities P R and P A . There are two bounds: (i) the minimum probability required to accept evidence that all m traits are regionally associated P R * and (ii) the minimum probability required to accept that the causal variantterm·key for all m traits aligns at a single variant P A * . The BB algorithm accepts evidence supporting all m traits sharing a single causal variantterm·key if P R P A ≥ P R * P A * , after which the algorithm returns the HyPrColocterm·mth estimate of PPFC and stops. If either P R < P R * or P A < P A * there is insufficient evidence supporting all traits sharing a causal variantterm·key and the BB algorithm moves to the branch with maximum evidence supporting m − 1 traits sharing a causal variantterm·key. At this point the traits are partitioned into two clusters: one containing m − 1 traits deemed most likely to share a causal variantterm·key and a second cluster containing the remaining trait. We repeat this process of branch selection and partitioning on the cluster of m − 1 traits until we identify either: (A) a cluster of traits of size k ≥ 2 whose regional and alignment statistics satisfy P R P A ≥ P R * P A * , or (B) there is one trait left in the cluster. In scenario A, the HyPrColocterm·mth posterior probability that all k traits colocalize is presented and the remaining m – k traits are assessed for evidence of colocalizationterm·key using the branch selection and partitioning scheme. In scenario B, the trait is deemed not colocalize with any other trait in the sample and the BB selection algorithm is repeated using m − 1 traits. The entire process is repeated until all clusters of colocalized traits, whereby each cluster of traits colocalize at a distinct causal variantterm·key, have been identified, all other traits are deemed not to share a causal variantterm·key with any other trait.

Simulation study To create genomic loci with realistic patterns of LDterm·dom, for each simulation scenario we simulated 1000 datasets and for each dataset we resampled phased haplotypes from the European samples in 1000 Genomes 15 and randomly chose one of the first 50 regions confirmed to be associated with CHDterm·dom 16 . After removing variants with low MAF, i.e. MAF < 0.05, the number of SNPs analysed in these regions ranged from 228, in the APOE region, to a maximum of 1918 SNPs in the PDGFD region. The mean number of SNPs was 881.6. Unless stated otherwise, for traits that have a causal variantterm·key in the region, the variant explains 1% of trait variance. To go some way to mirroring real analyses, each trait was assumed to be measured in studies with different sample sizes, i.e. the sample size for the i -th study ( N i ) was randomly chosen from the set N i ∈ {1,000, 5,000, 10,000, 15,000, 20,000}. Variant-level priors were chosen for the simulation study: we set p = 10 −4 as in refs. 2 , 14 and, to assess sensitivity of results to the choice of conditional colocalizationterm·key prior p c , we ran each simulation three times for each of p c ∈ {0.05, 0.02, 0.01}. Note that p c = 1 − γ , so this is equivalent to γ ∈ {0.95, 0.98, 0.99}. For a pair of traits, colocalizationterm·key between the traits is 5 times more likely a-priori when setting p c = 0.05 relative to p c = 0.01. In the analysis of ten traits, however, colocalizationterm·key between all ten traits is around 1 million times more likely a-priori when setting p c = 0.05 relative to p c = 0.01. The prior probability of colocalizationterm·key is still very small ~10 −11 when setting p c = 0.05, however. Hence, the different values of p c we have chosen can result in substantial differences in the prior probability of colocalizationterm·key.

Violations of the single causal variantterm·key assumption These data were generated under three scenarios: (i) a single cluster of 10 colocalized traits, each trait shares a single causal variantterm·key and 5 traits have secondary distinct causal variantsterm·key; (ii) a single cluster of 6 colocalized traits, each of the 6 traits share a single causal variantterm·key and 3 traits have secondary distinct causal variantsterm·key, the remaining 4 traits do not have causal variantsterm·key and; (iii) 2 clusters of colocalized traits, cluster 1 comprises 6 traits sharing a single causal variantterm·key with 3 of 6 traits having secondary distinct causal variantsterm·key, cluster 2 comprises 4 traits sharing a single causal variantterm·key with 2 of 4 traits having secondary distinct causal variantsterm·key. To maximize the number of traits with additional causal variantsterm·key in a cluster (up to the maximum of 5), in scenarios (ii) and (iii) the total number of clusters of colocalized traits were reduced relative to the single causal variantterm·key assessment.

Measuring the accuracyterm·mth, true positive and false positive ratesterm·mth of HyPrColocterm·mth Accuracyterm·mth = T P + T N T P + T N + F P + F N , True positive rateterm·mth T P R = T P T P + F N , False positive rateterm·mth ( F P R ) = F P F P + T N where TP and TN denote the true positive count and true negative count, and FP and FN denote the false positive count and false negative count. Hence, accuracyterm·mth is the proportion of traits that are correctly identified as colocalizing or not colocalizing. To compare HyPrColocterm·mth with pairwise methods, we compute the TP, FP, TN and FN rates by aggregating information across all pairs of traits in the sample. A TP is measured when a pair of observations are correctly deemed to colocalize, a FP is measured when a pair of traits are incorrectly identified as colocalizing, a FN is recorded when a pair of traits are wrongly deemed not to colocalize and a TN is recorded when a pair of traits are correctly identified as not colocalizing.

When thresholding the posterior probability of colocalizationterm·key, the TP, FP, TN and FN rates are computed after excluding traits which do not to colocalize with any other trait such that P A P R > 0.7. In the simulation study which allowed each trait a maximum of one causal variantterm·key in the region and with respect to scenarios (i), (ii) and (iii), when setting p c = 0.05 HyPrColocterm·mth identified clusters of colocalized traits with P R P A > 0.7 in approximately 70%, 93% and 99% of simulated datasets, when p c = 0.02 in approximately 65%, 91% and 98% datasets, reducing to around 60%, 86% and 97% of datasets when p c = 0.01. Pairwise COLOC identified pairs of colocalized traits with P R P A > 0.7 in over 96% of simulated datasets, across all three scenarios and specifications of p c . In the simulation study which allowed a maximum of two causal variantsterm·key per trait, these figures reduced: when setting p c = 0.05 HyPrColocterm·mth identified clusters of colocalized traits with P R P A > 0.7 in approximately 65%, 80% and 93% of simulated datasets, when p c = 0.02 in ~60%, 72% and 92% datasets, reducing to around 55%, 65% and 85% of datasets when p c = 0.01. Pairwise COLOC identified pairs of colocalized traits with P R P A > 0.7 in over 94% of simulated datasets, across all three scenarios and specifications of p c Application to CHDterm·dom and cardiovascular risk factors The GWASterm·mth results used in the assessment of colocalizationterm·key of CHDterm·dom with related traits were taken from large-scale analyses of CHDterm·dom 17 , blood pressure ( http://www.nealelab.is/uk-biobank ), adiposity measures ( http://www.nealelab.is/uk-biobank ), glycaemic traits 18 , renal function 19 , type II diabetes 20 , lipid measurements 21 , smoking 22 , rheumatoid arthritis 23 and educational attainment 24 (Table S1 ). All datasets had either been imputed to 1000 Genomes 15 prior to GWASterm·mth analyses or were imputed up to 1000 Genomes from the summary results using DIST 48 (INFO > 0.8). We performed colocalizationterm·key analyses in two steps. In step one, we assessed colocalizationterm·key of CHDterm·dom with the 14 risk-factors in pre-specified LDterm·dom blocks from across the genome 25 . We used a conservative variant-level prior structure with p = 1 × 10 −4 and γ = 0.98, i.e. 1 in 500,000 variants are expected to be causal for two traits, and set strong bounds for the regional and alignment probabilities, i.e. P R * = P A * = 0.8 so that the algorithm identified a cluster of colocalized traits only if P R P A > 0.64. The full results from this analysis are available at https://jrs95.shinyapps.io/hyprcolocterm·mth_chdterm·dom .

To prioritise candidate causal genes in regions where CHDterm·dom and at least one related trait colocalized, we re-ran the colocalizationterm·key analysis and included whole blood cis eQTL 27 (31,684 samples) and cis pQTL 28 (3301 samples) data in addition to the primary traits in a second step, using the same LDterm·dom blocks as before. A colocalizationterm·key analysis was performed for every transcript with data within each region. cis eQTL were defined 1MB upstream and downstream of the centre of the gene probe (1828 genes were analysed across the 43 regions). cis pQTL were defined 5MB upstream and downstream of the transcript start site (854 proteins were analysed across the 43 regions). We integrated gene expression information taken from whole blood tissue as: (i) the eQTLGen dataset 27 has a large sample size relative to other publicly available gene expression data resources and; (ii) the pQTL data were also measured in whole blood tissues, so there was consistency in the tissue analysed.

Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Supplementary information Supplementary Information Peer Review File Reporting Summary