Causal inference for heritable phenotypic risk factors using heterogeneous genetic instruments.

Over a decade of genome-wide association studies (GWAS) have led to the finding of extreme polygenicity of complex traits. The phenomenon that "all genes affect every complex trait" complicates Mendelian Randomization (MR) studies, where natural genetic variations are used as instruments to infer the causal effect of heritable risk factors. We reexamine the assumptions of existing MR methods and show how they need to be clarified to allow for pervasive horizontal pleiotropy and heterogeneous effect sizes. We propose a comprehensive framework GRAPPLE to analyze the causal effect of target risk factors with heterogeneous genetic instruments and identify possible pleiotropic patterns from data. By using GWAS summary statistics, GRAPPLE can efficiently use both strong and weak genetic instruments, detect the existence of multiple pleiotropic pathways, determine the causal direction and perform multivariable MR to adjust for confounding risk factors. With GRAPPLE, we analyze the effect of blood lipids, body mass index, and systolic blood pressure on 25 disease outcomes, gaining new information on their causal relationships and potential pleiotropic pathways involved.

Introduction

Understanding the pathogenic mechanism of common diseases is a fundamental goal in clinical research. As randomized controlled experiments are not always feasible, researchers are looking towards Mendelian Randomization (MR) as an alternative method for probing the causal mechanisms of common diseases [1]. MR uses inherited genetic variations as instrumental variables (IV) to interrogate the causal effect of heritable risk factor(s) on the disease of interest. The basic idea is that at these variant loci, the inherited alleles are randomly transmitted from the parents to their offsprings according to Mendel’s laws. Thus, the genotypes are independent from non-heritable confounding variables which may obfuscate causal estimation in parent-offspring studies. More generally, such independence also approximately holds for population data such as those collected in genome-wide association studies (GWAS) when individuals share the same ancestry [2]. With the accumulation of data from GWAS, there is increasing interest in MR approaches, especially in approaches that only rely on GWAS summary statistics that are publicly available [2, 3].

How well Mendelian Randomization works depends on how well the genetic variant loci used as instruments abide by the rules of IV. These rules dictate that, if the genetic locus has an effect on the disease outcome, it should be only through pathways mediated by the risk factor(s) of interest. This rule, termed exclusion restriction, is violated when there is horizontal pleiotropy, defined as the case where the genetic variant can influence the disease through pathways other than the given risk factor(s) [4]. There has been much recent attention on this issue [5-16] in MR, yet our understanding is far from complete. Current methods rely on different assumptions on the pattern of horizontal pleiotropy, often driven by statistical convenience rather than what geneticists have learned from real data. What assumptions on pleiotropy and genetic effects would be suitable? Would it be possible to learn the degree of pleiotropy from the data? Could we perform model diagnosis utilizing only GWAS summary statistics?

The pleiotropy issue that muddles Mendelian Randomization studies is, in a large part, due to the fact that complex traits are extremely polygenic [15, 17–24]. Accumulating evidence from GWAS studies indicates that many complex diseases may have an omnigenic architecture where all genes affect every complex trait [25]. While a few genes might be “core” genes, a large proportion of genes may have non-zero effects on diseases and their risk factors. Thus, in an MR study, many genetic instruments, if not all, may affect the disease through their effects on other unmeasured risk factors. In other words, in an MR analysis, not only would we expect horizontal pleiotropy to be a pervasive issue across all genetic variants, any disease or complex risk factor would also be associated with a large number of SNPs across the whole genome. Many existing MR methods rely on the assumption that pleiotropic effects sparsely involve only a few SNPs, which directly counters these recent insights. Methods that don’t assume sparsity often require the instrument strength independent of direct effect (InSIDE) assumption [6], which can be rather optimistic. Recently, a few new methods relaxed the InSIDE assumption to consider “correlated pleiotropy” through one pleiotropic pathway [11, 12, 15, 26]. However, when pleiotropic pathways exist, there would often be an issue in identifying the true causal effect of the risk factor, and most methods are restricted to allow for only one pleiotropy pathway. Armed with these assumptions, most existing methods also utilize only the few SNPs that have the strongest association with the risk factor as instruments, ignoring the SNPs that are weakly associated. In this work, we will show that weakly associated SNPs are also informative, and that a model combining weak and strong SNPs can increase the accuracy and stability of our estimations in some scenarios.

We propose a comprehensive statistical framework for causal effect estimation under the realistic assumption that pleiotropy may be pervasive across the genome. The framework, called GRAPPLE (Genome-wide mR Analysis under Pervasive PLEiotropy), facilitates interactive identification of multiple pleiotropic pathways and incorporates all SNPs associated with the risk factor at varying p-value thresholds into the analysis. GRAPPLE builds upon a previous statistical framework we developed called MR-RAPS [10] under the InSIDE assumption, but is much more comprehensive and flexible. GRAPPLE emphasizes the detection of multiple pleiotropic pathways when the InSIDE assumption is violated as well as the determination of the causal direction. GRAPPLE further addresses two common challenges: how to jointly estimate the effects with multiple risk factors to reduce correlated pleiotropy, and how to integrate cohorts with overlapping samples. The estimation accuracy of GRAPPLE is examined through validations involving both real studies and simulations.

We apply GRAPPLE to investigate the causal effects of 5 risk factors (three plasma lipid traits, body mass index, and systolic blood pressure) on 25 common diseases. Although there have been several causal effect screening studies [9, 15, 26] for these risk factors and diseases, the multi-modality analysis enabled by GRAPPLE brings forth new insights on the pleiotropic landscape of these diseases and, thus, an improved understanding of the causal risk factors. Specifically, we will reexamine the role of lipid traits on coronary artery disease and type-II diabetes, where the results from MR studies have been under heated debate [2, 27, 28]. The R package GRAPPLE is publicly available for installation at https://github.com/jingshuw/GRAPPLE.

Results

Model overview

From the causal model to GWAS summary statistics

Our framework starts with a set of structural equations that jointly specify the generative model on the disease Y that relies on K observed risk factors = (X1, ⋯, X) of interest, and the vector = (Z1, Z2, ⋯) containing all genetic information of an individual (Fig 1a). Here represents unknown non-heritable confounding factors and and E are random noise acting on X and Y respectively. The parameter of interest, , quantifies the causal effect of the vector of risk factors on Y. Mendel’s laws of inheritance suggest that the genotypes are randomized during conception and are generally independent of the environmental factors (). The function f(, , E) represents the causal effect of unmeasured risk factors on Y, which can be heritable (contributed by ) or non-heritable (contributed by ). The non-parametric functions f(⋅) and g(⋅) allow interactions among SNPs in and variables () in their causal effects on and Y. Under this model, there is horizontal pleiotropy for a SNP j if Z has nonzero association with f(, , E). This is the case, for example, when Z acts on Y through a pathway affecting unmeasured risk factors, or when Z is in linkage disequilibrium (LD) with such a locus.

Fig 1

Model overview.

a, The causal directed graph represented by structural equations (1). b, The existence of a pleiotropic pathway 2 (purple) can result in multiple modes of the profile likelihood. c, Multi-modality of the profile likelihood can reflect causal direction. d, The work-flow with GRAPPLE.

Now consider the case where only GWAS summary statistics, i.e. the estimated marginal associations between each SNP j and the risk factors/disease traits, are available and there are in total p SNPs selected. Let Γ be the true association between SNP j and Y, and γ be the vector of true marginal associations between SNP j and . Later, we will denote their estimated values from GWAS summary statistics as and . Then, as shown in Materials and Methods, the model (1) results in the linear relationship where for binary Y, the parameter in model (2) is a conservatively biased version of in model (1). This relationship holds even when the functions f(⋅) and g(⋅) in (1) are not linear. Here, α is the marginal association between Z and f(, , E), representing the unknown horizontal pleiotropy of SNP j.

One can immediately see that identifying is impossible without further assumptions regarding α. Early MR methods such as IVW [5] made the assumption that all instruments are valid satisfying α = 0. Other methods such as Weighted Median [7] or MR-PRSSO [9] assume that α is sparsely nonzero. However, the no or sparse pleiotropy assumption follows from statistical convenience rather than biological insights. As discussed in Introduction, horizontal pleiotropy is pervasive for most complex traits. One assumption that allows pervasive pleiotropy is to assume the InSIDE assumption [6] where , or alternatively, the random effect model [10, 16] where for most genetic instruments. Unfortunately, the InSIDE assumption can be easily violated if the pleiotropic effects of selected genetic variants are driven by shared pleiotropic mechanisms.

Some more recent MR methods such as LCV [26], MRMix [12], Contamination mixture [13] and CAUSE [15] have noticed this limitation of the InSIDE assumption and allow a subset of the genetic instruments to be associated with a common hidden pleiotropic pathway. For instance, using the above notation, both CAUSE and MRMix assumed that when for the SNPs that violate the inSIDE assumption, their pleiotropic effects satisfy (when K = 1) where aγ represents the correlated pleiotropic effects due to a confounding pathway and . This is a more realistic assumption than InSIDE, though there would then be an issue in distinguishing the true causal effect β from the pleiotropic direction β + a. Allowing for only one pleiotropic pathway also makes the model restrictive for real datasets.

Identify multiple pleiotropic pathways and the direction of causality

The key idea underlying GRAPPLE is that multiple pleiotropic pathways can be detected by using the shape of the profile likelihood function under no pleiotropy assumption. This allows us to probe the underlying causal mechanism, without explicit assumptions on pleiotropic patterns (Fig 1b). When K = 1, the GWAS summary statistics reduce to the scalar and , with their standard errors σ1 and σ2. From the central limit theorem, the joint distribution of approximately follows a multivariate normal distribution where θ is a shared sample correlation that can be estimated as (see Materials and methods).

When there is no horizontal pleiotropy in the p selected independent genetic instruments (α = 0 for j = 1, 2, ⋯, p), the robust profile likelihood [10] is given by, where ρ(⋅) is the Tukey’s Biweight loss, or any other robust loss functions. As described with more details in Materials and Methods, the profile likelihood is obtained by profiling out nuisance parameters γ1, ⋯, γ in the full likelihood from (3), which is further robustified by replacing the L2 loss with Tukey’s Biweight loss to increase the sensitivity of mode detection. Under no pleiotropy or InSIDE assumption, this function l(b) should have only one mode near the true causal effect b = β.

Now consider the case where a second genetic pathway (Pathway 2) also contributes substantially to the disease, and some instrument loci are also associated with Pathway 2 (Fig 1b). In this scenario, SNPs that are associated with X only through Pathway 2 can contribute to a second mode in the profile likelihood at location β + κ/δ, where κ and δ quantifies the causal effect of Pathway 2 on Y and its marginal association with X, respectively (Materials and Methods). Similarly, multiple pleiotropic pathways generally result in multiple modes of l(b). Thus, we can use multiple modes in a plot of l(b) to diagnose the presence of horizontal pleiotropic effects that are grouped by different pleiotropic pathways.

The existence of pleiotropic pathways complicates MR and makes the causal effects of the risk factors potentially unidentifiable. Specifically, when Pathway 2 exists, the GWAS summary statistics alone cannot provide information to distinguish β from β + κ/δ. Instead of making further untestable assumptions such as one pathway “dominates” the other, when multiple modes are detected, we suggest that whenever multiple modes are detected, the investigator should try to find biomarkers for each mode and collect more GWAS data to adjust for confounding risk factors. Specifically, GRAPPLE facilitates this by identifying marker SNPs of each mode, as well as the mapped genes and GWAS traits of each marker SNP (see Materials and methods). This allows researchers to use their expert knowledge to infer possible confounding risk factors that contribute to each mode. With GWAS summary statistics of these confounding traits, GRAPPLE can perform a multivariable MR analysis assuming the InSIDE assumption applies for the remaining horizontal pleiotropic effects (Materials and Methods).

The detection of multiple modes can be also used to determine the causal direction (Fig 1c). If the wrong causal direction is specified in model (1) and Y is a cause of X, the genetic variants associated with X can be classified in two groups: those associated with X through Y, and those associated with X through another pathway unrelated to Y. In the former case, γ = βΓ where β is the causal effect of Y on X, and these SNPs should contribute to a mode around 1/β. In the latter case, a SNP j satisfies γ ≠ 0 but Γ = 0, and would contribute to a mode of l(b) at 0. Thus, there will be two modes in the robust profile likelihood with one mode being around 0. This idea can be viewed as an extension of the bidirectional MR [29, 30]. Bidirectional MR is based on the assumptions that when X is a cause of Y, most of the genetic instruments for Y should be unassociated with X, because they affect Y through a different pathway, thus the reserve MR would indicate a zero effect of Y on X. GRAPPLE makes this inference more robust by making use of the fundamentally different shape of the robust profile likelihood plots in different directions. In the correct causal direction, the plot should only show one mode around the true causal effect β. In the incorrect reverse direction - when the true outcome is treated as the risk factor and the true risk factor is treated as the outcome - the plot of the robust profile likelihood will have two modes, one around 0, representing the variants directly related to the true outcome, and one around 1/β, representing the variants indirectly related to the true outcome through the true risk factor.

Weak genetic instruments: A curse or a blessing?

Besides the assumption of no-horizontal-pleiotropy, for a SNP to be a valid genetic instrument, it needs to have a non-zero association with the risk factor of interest. In most MR pipelines, SNPs are selected as instruments only when their p-values are below 10−8, which is required to guarantee a low family-wise error rate (FWER) for GWAS data. Using such a stringent threshold also helps to avoid weak instrument bias [31], where measurement errors in are not ignorable and lead to bias in . However, such a stringent selection threshold may result in very few, or even no instruments being selected with under-powered GWAS, and may still not be adequate to avoid weak instrument bias. Further, when our goal is to jointly model the effects of multiple risk factors (the setting where as a vector), it is unrealistic to assume that all selected SNPs have strong effects on every risk factor. In addition, the high polygenecity of complex traits indicates that the weak instruments far outnumbers strong instruments, and collectively, they may substantially improve the estimation accuracy.

In GRAPPLE, we use a flexible p-value threshold, which can be either as stringent as 10−8 or as relaxed as 10−2, for instrument selection. Based on the profile likelihood framework of MR-RAPS [10], GRAPPLE can provide valid inference for that avoids weak instrument bias for multiple risk factors even when the p-value threshold is as large as 10−2. This flexible p-value threshold is beneficial for several reasons. First, including moderate and weak instruments may increase power, especially for under-powered GWAS. Second, for MR with multiple risk factors where it is inevitable to include SNPs that have weak associations with some of the risk factors, we can obtain much more accurate causal effect estimations than methods that can only deal with jointly strong SNPs. More importantly, we can investigate the stability of the estimates across a series of p-value thresholds and get a more complete picture of the underlying horizontal pleiotropy. In practice, we suggest researchers to vary the selection p-value thresholds from a stringent one (say 10−8) to a relaxed one (say 10−2), both in the detection of multiple modes and in estimating causal effects.

The three-sample design to guard against instrument selection bias

The current two-sample design of MR uses one GWAS data for the risk factor and one for the disease. The selection of genetic instruments is performed with p-values reported in the GWAS data for the risk factor. However, selecting instruments from GWAS summary statistics can introduce bias, which is commonly referred to as the “winner’s curse”. Conditional on being selected, the magnitude of is generally larger than γ and introduces bias to the estimation of . When K = 1 where there is only one risk factor, the estimate will be biased towards 0, but there is no guarantee on the direction of the bias when K > 1. Among practitioners, a common belief is that the selection bias is negligible when only the strongly associated SNPs are selected as instruments.

However, this rule of thumb may not hold even when we only use that are genome-wide significant (p-value ≤10−8) (S1(a) Fig). Thus, we strongly advocate using a three-sample GWAS summary statistics design (Fig 1d). To avoid the selection bias, selection of genetic instruments is done on another GWAS dataset for the risk factor, whose cohort has no overlapping individuals with both the risk factor and disease cohorts. In addition, to simplify calculation and avoid bias due to different LD structure in heterogeneous populations, we use LD clumping [32] to select independent SNPs in GRAPPLE (see Materials and methods). The three-sample design will also avoid possible selection bias introduced during clumping.

Summarizing the above points, a complete diagram of the GRAPPLE workflow is shown in Fig 1d. A researcher may start with a single target risk factor of interest. The shape of the robust profile likelihood provides information on possible pleiotropic pathways. If only a single mode is detected, one can use GRAPPLE for the target risk factor. This is equivalent to using the original MR-RAPS. If multiple modes are detected, the researcher needs to seriously consider how to adjust for pleiotropic pathways. Researchers can use the marker SNP/gene/trait information that GRAPPLE provides to investigate each mode, decide on which confounding risk factors to adjust for, and collect extra GWAS data for them. GRAPPLE can then be used to jointly estimate the causal effects of the original and the additional risk factors.

Assessment of GRAPPLE with real studies

Combine weak and strong genetic instruments under no pleiotropy

We first examine whether GRAPPLE provides reliable statistical inference using instruments with different strength under an artificial setting with real GWAS summary statistics. In this setting, we make the “artificial risk factor” X and the “artificial disease” Y be the same trait from two non-overlapping cohorts, thus γ = Γ while for any SNP j. Though the structural equation describing the causal effect of X on Y is not well defined., the linear relationship model (2) from which we estimate β still holds with β = 1 and α = 0. In other words, we are not estimating a meaningful “casual” effect, but are in a special case where the true β is known. This setup can be used to verify the validity of MR methods under no pleiotropy.

Specifically, we consider three traits: Body mass index (BMI), Type II diabetes (T2D) and height from the GIANT and DIAGRAM consortia where sex-specific GWAS data are available [33-35]. The female cohort is used to get and the male cohort is used to get . As a three-sample design, the UK Biobank data for corresponding traits are used for SNP selection. If we assume that all selected instruments have no gender-specific association with the traits, the true β would equal 1. For benchmarking, we compare the performance of GRAPPLE with CAUSE [15] and three other widely adopted MR methods: inverse-variance weighted (IVW) [5], MR-Egger [6] and weighted median [7] with the same three-sample design.

We compare different p-value thresholds for instrument selection, ranging from a stringent threshold of 10−8 to a relaxed threshold of 10−2 (Fig 2a). GRAPPLE provides roughly unbiased estimates of β no matter which threshold is used, showing that it does not suffer from weak instrument bias. Surprisingly, the other MR methods are biased even with a stringent p-value threshold.

Fig 2

Performance evaluation.

a, Estimation of β across selection p-value thresholds under no pleiotropy. Error bars show 95% Confidence intervals and the numbers are the number of independent SNPs obtained at each threshold. b, Estimation of β across three non-overlapping categories of SNPs: “strong”, “moderate” and “weak”. The numbers are the number of SNPs in each category. c, Identifying causal directions by multi-modality with MR reversely performed. The selection p-value threshold is kept at 10−4. d, three modes detected in the profile likelihood with selection p-value threshold 10−4 for CRP on CAD. Marker genes and GWAS traits (in parenthesis) are shown for each mode. e, estimation of the CRP effect β at different p-value selection threshold with each method. The numbers are the estimated , with * indicating p-value below 0.05 and ** indicating p-value below 0.01.

Notice that for T2D, the confidence intervals of GRAPPLE do get narrower with increasing p-value thresholds (Fig 2a), showing the potential power gain of including weak instruments in less powerful GWAS studies. In addition, we simulate synthetic GWAS summary statistics of the risk factor and disease (see S1 Text for details) and confirm that the estimated β indeed gets more accurate with the inclusion of weakly associated SNPs (S2(a) Fig). In the simulations, we also use GRAPPLE to adjust for measured confounding risk factors and compare the performance with MVMR [36], a commonly used multivariable MR method. As discussed earlier, in multivariable MR, the inclusion of SNPs that are weakly associated with at least one risk factors is inevitable. As GRAPPLE does not suffer from weak instrument bias, we see that it provides accurate estimates of the causal effects as well as reliable confidence intervals with both stringent and mild p-value thresholds (S3–S5 Figs).

Finally, we demonstrate that to avoid bias, the three-sample design is necessary no matter which MR method is used. As shown in S1(a) Fig, the two-sample design where we use the same cohort of the risk factor for selection can result in biased casual effects estimation, and the bias occurs with most MR methods even when we only select the strongly associated SNPs.

Weak SNPs provide reliable causal estimates under pleiotropy

Next, we examine whether or not the weak instruments are more vulnerable to pleiotropy, which can be a concern for including the weak SNPs. We compare four risk factor and disease pairs that cover eight different complex traits, including the effect of BMI on T2D, low-density cholesterol concentrations (LDL-C) on coronary artery disease (CAD), height on smoking, and systolic blood pressure (SBP) on stroke (Fig 2b). The GWAS summary data are collected from the original study repositories [37-42].

We test whether independent sets of strongly and weakly associated SNPs can provide consistent estimates of the causal effects of the risk factors. SNPs passing the p-value threshold 10−2 in the cohort for selection are divided into three non-overlapping groups after LD clumping: “strong” (p ≤ 10−8), “moderate” (10−8 < p ≤ 10−5), and “weak” (10−5 < p ≤ 10−2). The SNPs across groups are used separately to obtain group-specific estimates of the causal effect β. We observe that for all the four pairs, the estimates are stable across groups (Fig 2b). Though the “weaker” SNPs provide estimates with more uncertainty due to limited power, the estimates are consistent with those from the “strong” group. Other MR methods also show some level of consistency in estimating β across different sets of instruments, but perform less well due to weak instrument bias (S1(b) Fig). To conclude, in the analysis of these four pairs of traits, we do not see any evidence that weakly associated SNPs provide more biased estimates than strong instruments due to horizontal pleiotropy. In contrast, as with the strong instruments, the weakly associated SNPs may also provide useful information to infer the causal effects of the risk factors.

Identify direction of causality for known causal relationships

We also examine the performance of GRAPPLE in identifying the causal direction with the shape of the profile likelihood. For the causal direction, we focus on the two pairs of traits with known causal relationship: BMI on T2D, and LDL-C on CAD. We switch the roles of the risk factor and disease to see if the correct direction can be revealed. Specifically, we treat T2D and CAD as the “risk factor”, and BMI and LDL-C as the corresponding “disease” (Fig 2c). For T2D, the cohort for the other gender is used for SNP selection and for CAD, the risk factor cohort used is from [43] and the selection p-values are from [44]. As expected, we see that when the roles of the risk factor and disease are reversed, the robust profile likelihood shows a main mode at 0, and a weaker mode around .

Detect multiple modes to identify pleiotropic pathways

Finally, we test the ability of GRAPPLE to identify multiple pleiotropic pathways with the analysis of the C-reactive protein (CRP) effect on CAD. C-reactive protein has been found to be strongly associated with the risk of heart disease while many SNPs that are associated with the C-reactive protein also seem to have pleiotropic effect on lipid traits [45]. Previous MR analyses only included SNPs that are near the CRP gene to guarantee a free-of-pleiotropy analysis [46, 47] and found that CRP has no causal effect on CAD. Now, instead of only using SNPs near the CRP gene, by using associated SNPs across the whole genome that are known to involve pleiotropy pathways, can GRAPPLE identify the existence of these pathways and still obtain the correct estimate of the C-reactive protein effect?

CRP GWAS data from [48] are used for selection and the data from [49] using a larger cohort is used for getting . Similar to a multi-modality pattern already reported in [11], our robust profile likelihood shows a pattern of three modes, indicating the existence of at least three different pathways (Fig 2d). One mode is negative, one is positive and the third is around zero. The negative mode involves a few marker genes including HNF1A and PVRL2, with a marker trait LDL-C. The positive mode has marker traits pulmonary function and the C-reactive protein, and the few markers genes (IL6R, ARHGAP10, BCL7B, PABPC4) are also involved in immune response and lung cancer progression [50, 51]. The mode at 0 has marker genes CRP and LEPR, and only one marker trait, C-reactive protein.

We compare across 3 p-value thresholds (10−8, 10−5, 10−3) and check how the existence of multiple pathways affects causal estimates of the effect of C-reactive protein in MR methods using SNPs across the genome. Including C-reactive protein as the only risk factor, all bench-marking methods give a negative estimate of the CRP effect, which is possibly driven by the bias from an LDL-C induced pleiotropic pathway (Fig 2e). MR-RAPS is the estimation method used in GRAPPLE if we only use one risk factor, and the three other bench-marking methods give incorrect inference of the CRP effect with a p-value of β below 0.01 for at least one SNP selection threshold (notice that the weak instrument bias is towards 0 as shown in Fig 2a, thus the significance at p-value threshold 10−3 for MR-Egger and IVW cannot be explained by weak instrument bias). In contrast, after using two risk factors: C-reactive protein and LDL-C, where LDL-C is an identified confounding risk factor from the marker SNPs in Fig 2d, the estimates of CRP effect are much closer to 0 compared with that without including LDL-C. This analysis illustrates how GRAPPLE can detect pleiotropic pathways, provide information to identify the confounding risk factors to adjust for, and obtain correct inference after adjusting for these risk factors.

As a complement to the analysis on CRP, we also use simulations (S1 Text) and generate synthetic disease traits to evaluate the precision and recall in the detection of multiple modes and marker SNPs when there are pleiotropic pathways. We consider scenarios with one or two pleiotropic pathways caused by hidden confounding risk factors, and vary the genetic correlations between these hidden factors and the target risk factor. A higher genetic correlation corresponds to a larger proportion of SNPs that have a correlated pleiotropic effect. We observe that the detection of multiple modes is most powerful when the genetic correlation is neither too large nor too small (S2(b) Fig). If the genetic correlation is too high, then there are not enough SNPs to contribute to the mode of the true causal effect, while if the genetic correlation is too low there will be too few SNPs to contribute to the pleiotropic modes. Including weaker SNPs will decrease the sensitivity in mode detection but can increase the recall of true marker SNPs. In our simulations, we also observe that all univariable MR methods can perform poorly in estimating the true causal effect in the presence of pleiotropic pathways (S3–S5 Figs).

A causal landscape from 5 risk factors to 25 common diseases

Finally, we apply GRAPPLE to interrogate the causal effects of 5 risk factors on 25 complex diseases. The five risk factors are three plasma lipid traits: LDL-C, high-density lipoprotein cholesterol (HDL-C), triglycerides (TG), BMI and SBP. The diseases include heart disease, Type II diabetes, kidney disease, common psychiatric disorders, inflammatory disease and cancer (Fig 3a). The GWAS summary statistics are from studies [35, 37, 38, 41, 42, 52–70] and downloaded from the GWAS catalog [71]. For each pair of the risk factor and disease, we compare across p-value thresholds from 10−8 to 10−2. As a summary of the results, Fig 3a illustrates the average number of modes detected across the p-value thresholds for SNP selection (for modes at each p-value threshold, see S6 Fig). Besides the number of modes, Fig 3a also shows the p-values for each risk factor when GRAPPLE is performed with only the single risk factor (see also S6 Fig and Materials and methods). These p-values are not valid when there are pleiotropic pathways.

Fig 3

Screening with GRAPPLE.

a, Landscape of pleiotropic pathways on 25 diseases. The colors show average number of modes across 7 different selection p-value thresholds. The “+” sign shows a positive estimated effect and “−” indicates a negative estimated effect, with the p-value for each cell a combined p-value (see Materials and methods) of replicability across 7 thresholds using the single risk factor. These p-values are not multiple-testing adjusted across pairs. b, Multi-modality of the profile likelihood for effect of HDL-C on CAD at 2 different selection p-value threshold. Vertical bars are positions of marker SNPs (), labeled by their mapped genes (only unique gene names are shown). c, Multivariable MR for the effect of 5 risk factors on CAD. d, Multivariable MR for the effect of 4 risk factors on CAD. The Error bars are 95% confidence intervals.

Fig 3a shows that multi-modality can be detected in many risk factor and disease pairs. Multi-modality is most easily seen using the stringent p-value threshold 10−8 (S6 Fig). However, we find that some modes are contributed by a single SNP thus is more likely an outlier than a pathway. For instance, the effect of stroke on LDL-C shows two modes when the p-value threshold is 10−8 or 10−7 (one mode around −2.3 and another mode near 0.08). However, the negative mode only has one marker SNP (rs3184504) which has been found strongly associated with hundreds of different traits according to GWAS Catalog while the other mode has hundreds or marker genes. After removing the SNP rs3184504, the mode disappears. Such a mode also disappears when we increase the p-value threshold to include more SNPs as instruments. Thus, the average number of modes serves as a strength of evidence for the existence of multiple pleiotropic pathways. When a risk factor and disease pair show multi-modality, the p-values from GRAPPLE using the single target risk factor are no longer valid, and the researchers need further investigations of the modes.

First, consider the well-studied, often-debated relationship between CAD and the lipid traits. All five risk factors show highly significant effects, though multi-modality is detected in HDL-C and SBP. In our results for HDL-C, with different p-value thresholds, three modes in total can show up, two being negative and one positive, indicating that the pathways from HDL-C to CAD is complicated (Fig 3b). Fig 3b shows that one negative mode is contributed by SNPs near genes LPL and BUD13, which are strongly associated with triglycerides. Another positive mode is contributed by SNPs near genes ALDH1A2 and PSKH1, which is related to respiratory diseases [72]. The markers of the other negative mode are mapped to genes including LIPG and CETP.

Since the effects of the lipid traits are generally complicated, we combine all 5 risk factors and run an MR jointly with GRAPPLE (Fig 3c) with different p-value thresholds. After adjusting for other risk factors, the two most prominent risk factors for heart disease are LDL-C and SBP, while the protective effect of HDL-C stays negligible, as does the risk conferred by TG. So these results show that HDL-C as a single measurement does not seem to have a protective effect on heart disease with multiple complex pathways linking HDL-C and heart disease. Researchers have suggested analyzing different subgroups of HDL-C as smaller particles tend to have a stronger protective effect [73].

Lipids are also involved in a number of biological functions including energy storage, signaling, and acting as structural components of cell membranes and have been reported to be associated with various diseases [74-77]. Besides CAD, another disease that most likely involves the lipid traits is the Type II diabetes (Fig 3a). T2D is associated with dyslipidemia (i.e., higher concentrations of TG and LDL-C, and lower concentrations of HDL-C), though the causal relationship is still unclear [78]. In the meantime, evidence has emerged that LDL-C reduction with statin therapy results in a modest increase in risk of T2D [74]. For the MR analyzing each risk factor alone, we see potential protective effects of LDL-C and HDL-C on T2D but also multi-modality patterns. Two modes show up in the profile likelihood from HDL-C to T2D where one negative mode has a marker gene LPL and a mode near 0 with marker genes CETP and AC012181.1. Thus, we include all 3 lipid traits, along with BMI and run a joint model for these 4 risk factors using GRAPPLE (Fig 3d). Our result indicates a mild protective effect of HDL-C and LDL-C on T2D, and close to the null but imprecise estimate for TG.

Discussion

We propose a comprehensive framework, GRAPPLE, that utilizes both strongly and weakly associated SNPs to understand the causal relationship between complex traits. GRAPPLE is robust to pervasive pleiotropy and can identify multiple pleiotropic pathways. The multivariable MR performed by GRAPPLE can adjust for known confounding risk factors.

GRAPPLE incorporates several improvements over existing MR methods. It avoids weak instrument bias by dealing with measurement errors of the SNP associations on the risk factors with profile likelihood. Our likelihood from (9) is similar to the likelihood used in [79], but our likelihood allows modelling pervasive pleiotropy as long as the InSIDE assumption holds for most SNPs. The multi-modality visualization shares similarities with [8], which estimates the causal effect by the global mode, but we provide a more comprehensive analysis to identify multiple pleiotropic pathways by the local modes. Our causality direction identification is related to bidirectional MR where they used the assumption that if we reverse the role of risk factor and disease, the estimated causal effect is likely to be 0. We use this idea in a more principled way and can avoid bias when SNPs affecting the disease through the target risk factors are also selected as variants for the disease in the reverse MR. Finally, as the intercept term in MR-Egger is not invariant to the arbitrary assignment of effect alleles for each SNP, indicating a deficiency of the method, GRAPPLE does not include any intercept term.

GRAPPLE needs a separate GWAS cohort of the exposure for SNP selection, which is necessary for valid inference with weakly associated SNPs. Actually, as shown in S1(a) Fig, the three-sample design is needed for other MR methods as well to avoid selection bias. In some domains, it is hard to obtain multiple good-quality public GWAS summary statistics with non-overlapping cohorts. We call for the release of stage-specific or study-specific GWAS data summary statistics to the public in the future.

In GRAPPLE, we still require using a p-value threshold, though it can be as relaxed as 10−2, instead of requiring no p-value threshold at all. There are two main reasons for this requirement. One consideration is to increase power, as including too many SNPs with γ = 0 or extremely small would instead increase the variance of [10, 80]. Another consideration is that we would not want unmeasured risk factors that are unassociated (or very weakly associated) with target risk factors to bring in large pleiotropic effects with SNPs that mainly affect these unmeasured risk factors. The chance of including these SNPs would be much lower by requiring a relaxed p-value threshold.

Finally, when discussing the causal effect of a risk factor, one implicit assumption we use is consistency, assuming that there is a clear and only one version of intervention that can be done on the risk factor. However, interventions on risk factors such as BMI are typically vague [81]. For instance, there can be multiple ways to change weight, such as taking exercise, switching to different diet or conducting a surgery. It is common sense that these different interventions would have different effects on diseases, though they may change BMI by the same amount. However, the basic MR principle of gene-environment equivalence [82] suggests that whilst there will be genetic mimics of increased physical activity and decreased calorie intake, there will be no such mimics for having surgery. Basic biological principles indicate which inferences can be sensibly made. For example, cholesterol has multiple functions in our bodies and is involved in multiple biological processes. Intervening in different biological processes to change the concentration of lipid traits may, in principle, have different effects on disease. However, for many LDL-C lowering drugs the direct genetic mimics produce effects as predicted by RCTs of these pharmaceutical agents [83], demonstrating that gene-environment equivalence applies. We think that our causal inference using GRAPPLE, along with the markers we detect, would provide abundant information to deepen our understanding of the risk factors. However, one still needs to be careful when giving causal interpretations of the results. One recommendation in practice is to triangulate the results from MR with other sources of evidence [84, 85].

Materials and methods

Model details

The structural Eq (1) where = (X1, X2, ⋯, X) and = (β1, β2, ⋯, β) describe how individual level data are generated. To link it with the GWAS summary statistics data, denote which is the true marginal association between a SNP Z and risk factor X and which is the marginal association between Z and the causal effects of unmeasured risk factors on Y, i.e. the horizontal pleiotropic effect of Z on Y given . Then we can rewrite the structural equations into the following linear models: where corr(Z, ϵ) = 0 for any k and is guaranteed by the definitions of γ and α. By replacing in (6) with (5), we get where and . As Corr(Z, e) = 0, we conclude that Γ also satisfies that

Thus, parameters Γ also represent true marginal associations between SNP Z and the disease trait. This is how we result in working with Eq (2).

When the disease is a binary trait, the structural equation of Y changes to With the same argument, we have

If we further assume that for each genetic instrument j, Z is actually independent of e (instead of just being uncorrelated), then the odds ratio that is estimated from the marginal logistic regression will be approximately Γ/c with a constant c > 1 determined by the distribution of e. In other words, for binary disease outcomes, Eq (2) is still approximately correct with the in (2) being a conservatively biased (by a ratio of 1/c) version of the in (7) (for a detailed calculation, see A.1 of [10]).

GWAS summary statistics from overlapping cohorts

The GWAS estimated effect sizes (log odds ratios for binary traits) of SNP j are for the disease and a length K vector for the risk factors. As shown in [86] and derived in S2 Text, for any risk factor k we have where N and N are the total sample sizes for the disease and kth risk factor. N is the number of shared samples. The correlation of X and Y of any shared sample is Corr[Y, X]. Eq (8) shows that all the SNPs share the same correlation. As a consequence, we assume where Σ is the unknown shared correlation matrix.

Estimate the shared correlation Σ

To estimate Σ from summary statistics, we can use Eq (8). We first need to choose SNPs where γ = 0 for all risk factors k so that we can estimate the shared correlation using the sample correlation of the chosen SNPs. We choose all SNPs whose selection p-values p ≥ 0.5 for all k.

For these selected SNPs, denote the Z-values of for j = 1, ⋯, T as matrix Z where T is the number of selected SNPs. Then Σ is estimated as the correlation matrix of Z.

Instruments selection using LD clumping

In GRAPPLE, we need to first select a set of SNPs as genetic instruments to estimate the causal effects . Here, we only select independent SNPs to simplify the calculation. Besides the independence requirement, we only include SNPs that pass a p-value threshold to reduce the inclusion of false positives that can decrease power. To avoid selection bias, a separate cohort for each risk factor is used where the reported p-values in that cohort are used for instruments selection. Denote the selection p-value for SNP j and risk factor k as p, for multiple risk factors and a given selection threshold, we require the Bonferroni combined p-values K min(p) to pass the threshold. After that, we use LD clumping with PLINK [87] to select independent genetic instruments. The LD r2 threshold for PLINK is set to 0.001.

Estimate the effects β

Here, we perform statistical analysis assuming α ∼ N(0, τ2) for the pleiotropic effects, while robust to outliers where the pleiotropic effects for a few instruments are large.

Under model (9), Eq (2) and given Σ, the log-likelihood with GWAS summary statistics satisfy: up to some additive constant. Here, e = (1, 0, ⋯, 0).

Define for each SNP j the statistics where is the variance of and is the covariance between and in Σ. Then the profile log-likelihood that profile out parameters results in

As discussed in [10], maximizing would not give consistent estimate of τ2. Because of this and the goal of making robust to outlier SNPs with large pleiotropic effects, our optimization function is the adjusted robust profile likelihood defined as where ρ(⋅) is some robust loss function. By default, GRAPPLE uses the Tukey’s Biweight loss function: where c is set to its common default value 4.6851. We maximize (11) with respect to as well as solving the following estimating equation for the heterogeneity τ2 which is where with . The estimating equation satisfies at the true values of and τ2, thus can result in consistent estimate of τ2. For the details of estimating and τ2 as well building confidence intervals for them, see S2 Text.

Identify pleiotropic pathways via the multi-modality diagnosis

We use the mode detection of the robust profile likelihood (11) to detect multiple pleiotropic pathways. To increase sensitivity, we set τ2 = 0 and reduce the tuning parameter in the Tukey’s Biweight loss function to c = 3. Here we present a detailed argument on why mode detection can identify pleiotropic pathways.

If there is a confounding Genetic Pathway 2 , as shown in Fig 1a, that are missed, then we have the structural equation and also the linear model for a SNP j that only associates with Genetic Pathway 2 and uncorrelates with X conditional on . Similar to (5), we have

Plug in (13), we have

Thus, if there are enough SNPs like SNP j, they would contribute to another mode of (4) at β + κ/δ.

The same argument works for identification of the causal direction. Say there is another that affects Y but is uncorrelated with the risk factor X (δ = 0). The existence of such is common, unless X is the only heritable risk factor of Y. SNPs strongly associated with would not likely be selected when X is the exposure while would appear when the roles of X and Y are switched. These SNPs can be used to identify the causal direction, as as in the reverse MR, they contribute to a mode at 0, while the SNPs that affect Y through X will contribute to a mode at 1/β.

Select marker SNPs and genes for each mode

GRAPPLE uses LD clumping with a stringent r2 (= 0.001) threshold to guarantee independence among the genetic instruments. However, marker SNPs are not restricted to these independent instruments in order to get more biological meaningful markers. Marker SNPs are selected from a SNP set where the SNPs are selected using LD clumping with r2 threshold 0.05.

Assume that there are M modes detected at positions β1, β2, ⋯, β. Define the residual of SNP j () for mode m as where t(⋅, ⋅) is defined in Eq (10). SNP j is selected as a marker for mode m if |r|>t1 for any m′ ≠ m and |r|≤t0. By default, t1 is set to 2 and t0 is set to 1 which gives reasonable results in practice. When the marker SNPs are selected, GRAPPLE further map the SNPs to ENCODE genes where the marker SNPs locate and search for the traits that these SNPs are strongly associated with in GWA studies by querying HaploReg v4.1 [88] using the R package HaploR. The ratios of the marker SNPs are also returned for reference (shown as the vertical bars in Fig 3b).

Compute replicability p-values across SNP selection thresholds

Each p-value shown in Fig 3a summarizes a vector of p-values across 7 different selection p-value thresholds ranging from 10−8 ot 10−2 for each risk factor and disease pair. It reflects how consistent the significance is across SNP selection thresholds. Specifically, it is the partial conjunction p-value [89] for rejecting the null that β is non-zero for at most 2 of the selection thresholds. For a risk factor and disease pair k, let the p-values computed by using SNPs selected with the 7 thresholds p where s = 1, 2, ⋯, 7. Then rank them as p ≤ p ≤ ⋯ ≤ p, the partial conjunction p-value for the pair k is computed as 5p.

Additional evaluation results with real data.

a, Selection bias in MR methods when SNP selection and are obtained from the same GWAS dataset. True β ≈ 1 and error bars show 95% confidence intervals. The numbers are the number of clumped SNPs at different threshold. b, The estimate of β across three independent categories of SNPs with different association strengths for four risk factor and disease pairs using three other bench-marking MR methods. The numbers are the number of SNPs in each category, separated by the values of their selection p-values (dashed vertical lines).

(TIFF)

Click here for additional data file.

Simulation results.

a, Boxplots of the estimated β1 using different MR methods over 100 repeated random experiments when there are no correlated pleiotropy. We compare across three different β1 values (0.2, 0.5 and 1) with SNPs selected by three different selection thresholds: 10−8 for the top 106 SNPs, 10−5 for the top 217 SNPs and 0.01 for the top 422 SNPs. b, Performance of GRAPPLE in detecting multi-modality. In each setting with pleiotropic pathways, we evaluate three metrics: the detection rate of multi-modality, the precision of the identified marker genes of the pleiotropic pathways and the recall of true marker genes that are identified. Each color represent a different metric and each shape is for a different selection threshold. The title of each plot shows (β1, ⋯, β) in each setting where β1 is the true causal effect, and in each setting, we vary the genetic correlation between each genetic confounding risk factor and the risk factor of interest.

(TIFF)

Click here for additional data file.

Comparison of different MR methods in settings with pleiotropic pathways when the selection threshold is 10−8 (top 106 SNPs).

a, Boxplots of the estimated β1 using different MR methods over 100 repeated random experiments. b The actual coverage of the 95% confidence intervals of β1 provided by different methods. For CAUSE, we report the coverage of the 95% credible intervals of β1. The red dotted line shows the expected 0.95 nominal level. For the three settings in the second row with β1 = 0, the CI coverage is the same as 1−type I error.

(TIFF)

Click here for additional data file.

Same as S3 Fig with the selection threshold being 10−5 (top 217 SNPs).

(TIFF)

Click here for additional data file.

Same as S3 Fig with the selection threshold being 10−2 (top 422 SNPs).

(TIFF)

Click here for additional data file.

Additional results on identifying pleiotropic pathways on 25 diseases.

Each figure is for results obtained using one of the 7 p-value thresholds. The colors show the number of detected modes. The “+” sign shows a positive estimated effect and “−” sign shows a negative estimated effect.

(TIFF)

Click here for additional data file.

Simulation details.

(PDF)

Click here for additional data file.

Additional mathematical details of the statistical analysis.

(PDF)

Click here for additional data file.

A list of resources for GWAS datasets used in the paper.

(PDF)

Click here for additional data file.