Exploring causality with biliary atresia at different levels: two-sample Mendelian randomization study

Research design

Based on the two-sample MR analysis, we assessed 41 inflammatory cytokines, 731 immune cell traits, and 1400 metabolites for their potential causal relationships with BA as the outcome. The instrumental variables (IVs) used in MR analysis must adhere to three assumptions: (1) IVs must be exposure related, (2) IVs must be independent of other confounders, (3) IVs are unrelated to the outcome and can influence the outcome only through exposure.6 Figure 1 provides an overview of this MR study, and the study methods align with the Strengthening the Reporting of Observational Studies in Epidemiology-MR checklist.7

Request permission

GWAS data sources

We identified genetic predictors of 41 inflammatory cytokines from a comprehensive cytokine-related GWAS meta-analysis involving three independent cohorts. The study encompassed 8293 Finnish participants from the Cardiovascular Risk in Young Finns Study and the ‘FINRISK’ Studies (FINRISK1997 and FINRISK2002).8

Each immune trait from the GWAS data can be accessed publicly via the GWAS Catalog (https://www.ebi.ac.uk/gwas/, GCST0001391-GCST0002121).9 A total of 731 immune cell traits, comprising median fluorescence intensities reflecting surface antigen levels (n=389), relative cell (RC) counts (n=192), absolute cell (AC) counts (n=118), and morphological parameters (n=32), were included (online supplemental table 1). The GWAS for this immunologic profile used data from 3757 European individuals with no overlapping cohorts.

The GWAS dataset for 1400 metabolites was sourced from the study conducted by Chen et al.10 This represented the most comprehensive analysis of human metabolites, with complete summary statistics available through the GWAS Catalog (https://www.ebi.ac.uk/gwas/, GCST90199621-GCST90201020) in the public domain. The extensive GWAS involved 1091 metabolites and 309 metabolite ratios in 8299 individuals from the Canadian Longitudinal Study on Aging cohort. The dataset comprises 1091 metabolites, including 850 known metabolites and 241 unknown metabolites (online supplemental table 2). These known metabolites, categorized into peptide, nucleotide, amino acid, carbohydrates, cofactors and vitamin, energy, lipid, and xenobiotics metabolism, are based on the Kyoto Encyclopedia of Genes and Genomes database.

The GWAS dataset for BA was obtained from the study conducted by Chen et al.11 The first GWAS was conducted in a European-American cohort comprising 343 patients with isolated BA and 1716 controls. The second GWAS involved a European-American cohort with 156 patients presenting with BA and other extrahepatic abnormalities, alongside 212 controls. GWAS summary data can be accessed at https://www.ebi.ac.uk/biostudies/ (S-BSST182).

Selection of IVs

The IVs selected for this MR analysis are based on three fundamental assumptions. Initially, single nucleotide polymorphisms (SNPs) with association thresholds at p<1×10⁻⁵ were extracted for each exposure. Subsequently, the clumping program within PLINK software (V.1.90) was employed to prune these SNPs, using a linkage disequilibrium (LD) r² threshold of <0.1 within a 500 kb distance. The LD r² was calculated based on the 1000 Genomes Project as a reference panel.12 Third, to validate the second hypothesis, we employed the ‘phenoscanner’ function to control for confounding factors. Fourth, to ensure the validity of the third hypothesis, we excluded SNPs strongly associated with the outcome variable (p <1×10^–5). Additionally, alleles with palindromic effects were removed, and the explained variance (R²) and F-statistics parameters were used to assess the strength of IVs. A threshold of F-statistics >10 is generally recommended for MR analysis. Finally, we conducted MR analysis on metabolites comprising more than two SNPs, calculating R² and F-statistics as follows:

Where β represents the effect size of the genetic variant of interest; EAF denotes the effect allele frequency of the genetic variant of interest; se (β) indicates the SE of the effect size of the genetic variant of interest; N stands for the exposure sample size; and k is the number of SNPs.

MR analysis

This MR analysis primarily estimated the causal association using standard inverse variance weighting (IVW) methods.13 When the IVs satisfy all three assumptions of the valid IVs, the IVW approach can provide a consistent assessment of the causality of the exposure. However, if some IVs deviate from these assumptions, the analysis may yield inaccurate results. Hence, the following analyses were conducted: (1) heterogeneity among selected IVs was assessed using the Cochran Q-test with the corresponding p value. If the null hypothesis was rejected, a random-effect IVW was applied instead of a fixed-effect IVW.14 (2) Horizontal pleiotropy was evaluated through the intercept using the MR-Egger method to ensure that genetic variants were independently associated with both exposure and outcome.15 (3) Additional analyses using MR methods with different modeling assumptions and strengths (MR-Egger, simple mode, weighted median16, and weighted mode17) were performed to increase the stability and robustness of the results. (4) The MR pleiotropy residual sum and outlier (MR-PRESSO) method was used to identify and exclude horizontal pleiotropic outliers that might significantly influence the estimated results.18 (5) A ‘leave-one-out’ sensitivity analysis was employed to assess whether the results were influenced by any single SNP.15 Additionally, we conducted MR-Steiger and reverse MR tests to ascertain whether our results supported our hypothesis.

Evaluation of genetic correlation and directionality

Although SNPs associated with BA were excluded in the selection of IVs, unrelated SNPs might have mediated the inheritance of BA. Therefore, to ensure that causal effects were not confounded by exposure consistency with the outcome, the SNP-based Linkage Disequilibrium Score (LDSC) was used to estimate coinheritance by Χ² statistics for both traits.

Confounding analysis

Although we evaluated horizontal pleiotropy in the MR results through various sensitivity analyses to detect SNPs violating the MR assumption, residual confounding SNPs might still exist. The ‘phenoscanner’ function was employed to assess whether each SNP was associated with known risk factors for BA, such as viral infection and maternal exposure to certain chemicals.19 If any SNP was found to be associated with the mentioned confounders (p<1×10^–5), the MR analysis was re-executed after excluding these SNPs to validate the reliability of the results.

Statistical analysis

Statistical analyses were conducted using R V.4.1.2 software (http://www.Rproject.org). MR analysis was carried out using the ‘TwoSampleMR’ package (V.0.5.7).20 IVW, weighted median-based, and model-based approaches were performed with the use of the ‘MendelianRandomization’ package (V.0.7.0).21 The false discovery rate (FDR)-adjusted p value, below 0.2, suggests statistical significance.22

Causal assessment of 41 inflammatory cytokines for BA

To investigate the causal relationship between inflammatory cytokines and BA, a two-sample MR analysis was conducted. Using IVW as the primary method, three suggestive associations were detected following FDR adjustment (p_FDR<0.2); we identified three inflammatory cytokines as risk factors for the development of BA (figure 2A): eotaxin (odds ratio (OR)=1.45, 95% confidence interval (CI): 1.08 to 1.95, p_FDR=0.18, p=0.014), G-CSF (OR =4.21, 95% CI: 1.75 to 10.13, p_FDR=0.05, p=0.001) and MCP-1/MCAF (OR=1.53, 95% CI: 1.12 to 2.10, p_FDR=0.14, p=0.007).

Request permission

Causal assessment of 731 immune cell traits for BA

To investigate the causal relationship between immune cell traits and BA, a two-sample MR analysis was employed. Using IVW as the primary analysis method, a total of 34 suggestive associations were identified, following FDR adjustment (p_FDR<0.2); we identified three immune cell traits as risk factors for the development of BA, two of which were derived from RC counts and one from AC counts (figure 2B): CD8dim NKT/T cells ratio (OR =0.59, 95% CI: 0.45 to 0.77, p_FDR=0.06, p=0.000), CD8dim NKT counts (OR =0.58, 95% CI: 0.43 to 0.78, p_FDR=0.06, p=0.000), CD8dim NKT/lymphocyte ratio (OR=0.63, 95% CI: 0.49 to 0.81, p_FDR=0.06, p=0.000).

Causal assessment of 1400 metabolites for BA

To investigate the causal relationship between metabolites and BA, a two-sample MR analysis was conducted. Employing IVW as the primary analysis method, a total of 78 suggestive associations were identified following FDR adjustment (p_FDR<0.2); we identified only one metabolite as a risk factor for BA development, and it came from the category of unknown metabolites (figure 2C): X-12261 levels (OR=2.86, 95% CI: 1.73 to 4.74, p_FDR=0.06, p=0.000).

The detailed data of IVs were shown in online supplemental table 3. Generally, when statistical significance was observed in two other MR tests (p<0.05), the causal association is deemed robust. Notably, the association of MCP-1/MCAF, CD8dim NKT counts, CD8dim NKT/T cells ratio, CD8dim NKT/lymphocyte ratio, and X-12261 levels was deemed robust (online supplemental table 4). Although eotaxin and G-CSF were significant only in the IVW and weighted median methods (p<0.05), they still signified potential causal associations (online supplemental table 4).

Heterogeneity test

Heterogeneity was tested using the Cochran Q-test. The p values of seven traits (eotaxin, G-CSF, MCP-1/MCAF, CD8dim NKT count, CD8dim NKT/T cell ratio, CD8dim NKT/lymphocyte ratio, and X-12261 level) were greater than 0.05, indicating that no significant heterogeneity was found (online supplemental table 4).

Sensitivity analysis

Sensitivity analyses were performed to identify horizontal pleiotropy and potential outliers. All associations were scrutinized for potential horizontal pleiotropy using the MR-Egger intercept and MR-PRESSO global test. All p values in the results were greater than 0.05, indicating that there was no significant horizontal pleiotropy (online supplemental table 4). Additionally, we visualized scatter plots (figure 3), funnel plots (online supplemental figure 1), and ‘leave-one-out’ plot (online supplemental figure 2) from which no potential outliers and the possibility of horizontal pleiotropy were found.

Request permission

Reverse causality analysis

To test hypothesis 3 (‘IVs are unrelated to the outcome and can only influence the outcome through exposure’), MR-Steiger and reverse MR tests were conducted. Results indicated that in reverse causality analysis with BA as the exposure factor, the p values of the MR-Steiger for eotaxin, G-CSF, MCP-1/MCAF, and X-12261 levels were all below 0.05. And their p values of the reverse MR test were all above 0.05, confirming the study’s correct directionality (online supplemental tables 4 and 5). However, due to the lack of data for immune cell traits in the reverse MR test, their association with BA was verified solely through the MR-Steiger method. CD8dim NKT counts, CD8dim NKT/T cells ratio, and CD8dim NKT/lymphocyte ratio all had the p values of the MR-Steiger below 0.05, validating the study’s directionality (online supplemental table 4).

Linkage Disequilibrium Score

The results of LDSC showed that eotaxin (R_g=0.068, SE=0.110, p=0.537), G-CSF (R_g=0.572, SE=0.313, p=0.068), MCP-1/MCAF (R_g=−0.077, SE=0.081, p=0.342), and X-12261 levels (R_g=−0.060, SE=0.161, p=0.711), which indicated that the MR estimates were not disturbed by the shared genetic component (table 1). We estimated the SNP–SNP heritability (proportion of variance attributed to genome-wide SNPs) for the four associations using LDSC. The SNP heritability of metabolites varied from 0.015 (G-CSF) to 0.167 (MCP-1/MCAF) (online supplemental table 6). The immune cell traits dataset was not analyzed as it lacked the necessary data for LDSC completion.

Confounding analysis

Although sensitivity analysis in this study was conducted to assess SNP estimates, to meet assumption 2 (‘IVs must be independent of other confounders’), we examined whether all SNPs associated with the seven traits were independent of BA risk factors (viral infection and maternal exposure to certain chemicals). Results indicated that none of the seven traits were associated with any confounding factors (online supplemental table 7).