The causal relationship between obstructive sleep apnea and rheumatic disease: A bidirectional Mendelian randomization study

Study Design

We conducted a two-sample MR analysis to explore the causal relationship between OSA and rheumatic disease, including AS, RA, SLE, gout, SS, BD, PsA, and SSc. Genetic variation as an instrumental variable (IV) needs to follow three assumptions: (1) reliably and robustly correlated with exposure, (2) independent of risk factors and confounding factors, and (3) genetic variation is only associated with outcome through selective exposure rather than through alternative pathways. This analysis utilized recently published genome-wide association studies (GWASs) based on the summary data. This study aims to explore the causal relationship between OSA and rheumatic disease through univariate bidirectional two-sample MR analysis. We evaluated whether a causal relationship exists between OSA and rheumatic disease and investigated whether there is a causal relationship between rheumatoid immune diseases and OSA. The design of our MR framework is illustrated in Figure 1.

Figure 1

Flowchart of the causal inference between OSA and rheumatic disease. GWAS, genome-wide association studies; LD, linkage disequilibrium; SNPs, single nucleotide polymorphisms; MR, mendelian randomized; IVW, inverse variance weighted; WM, weighted median; OSA, obstructive sleep apnea; AS, ankylosing spondylitis; RA, rheumatoid arthritis; SLE, Systemic lupus erythematosus; SS, Sjogren’s syndrome; BD, Behcet disease; PsA, Psoriatic arthritis; SSc, Systemic sclerosis.

Data Sources

Summary statistics of OSA in the European population were obtained from FinnGen Release 9 (https://www.finngen.fi/en). These GWAS data included 38,998 OSA cases and 336,659 normal controls of European ancestry (375,657 subjects in total). Data for AS, gout, SS, BD, PsA, and SSc were extracted from the FinnGen Release 9 (https://www.finngen.fi/en). The data for SLE and RA were acquired from the Integrative Epidemiology Unit (IEU) Open GWAS database (https://gwas.mrcieu.ac.uk/). The sample sizes for each of rheumatic diseases were as follows: AS (2,860 cases and 270,964 controls), RA (14,361 cases and 42, 923 controls), SLE (5,201 cases and 9,066 controls), gout (8,489 cases and 240,862 controls), SS (2,495 cases and 365,533 controls), BD (85 cases and 365,522 controls), PsA (3,186 cases and 240,862 controls), and SSc (619 cases and 365,513 controls). Detailed information on genetic data is provided in Table 1.

Selection of Instrumental Variables (IVs)

In the current study, we classified single nucleotide polymorphisms (SNPs) as IVs. To fulfill three assumptions for the MR study design, firstly, we identified candidate IVs using summary statistics from GWAS with a significance threshold of P < 5e^-8 for exposure-associated SNPs. In order to obtain a relatively large number of IVs for BD, we selected SNPs on a lower significant threshold (P < 5e^-5). For SSc and PsA, we selected SNPs on a lower significant threshold (P < 5e^-6). We applied a clumping technique, filtering SNPs with a linkage disequilibrium (LD) threshold of R² < 0.001 and a window size of 10,000 kb to ensure independence among the IVs. Furthermore, to assess the strength of the selected IVs, we calculated the R² and F statistics for each SNP using the equation: R² = 2 × EAF × (1-EAF) ×β² and F statistic= R²× (N-2)/(1-R²), where EAF represents the effect allele frequency, β is an estimate of the genetic effect of each SNP on the exposure, and N is the sample size of the exposure GWAS data. An F value greater than 10 indicates that the SNP is sufficiently robust to mitigate potential bias.^[22] In addition, SNPs closely associated with OSA were eliminated by searching PhenoScanner (http://www.Phenoscanner.cam.ac.uk/) and GWAS Catalog (P < 1e^-5), such as obesity, fluid retention, adenotonsillar hypertrophy, smoking, and drinking. In the same way, we eliminated SNPs that may be associated with rheumatic disease, such as obesity, body mass index, smoking, and drinking. We carefully screened and removed these SNPs to ensure the integrity of the IVs. Following rigorous selection criteria, the remaining SNPs were considered eligible IVs.

Based on these rigorous criteria, 12 independent SNPs for OSA, 12 SNPs for AS, 55 SNPs for RA, 40 SNPs for SLE, 8 SNPs for gout, 6 SNPs for SS, 17 SNPs for BD, 7 SNPs for PsA and 9 SNPs for SSc were separately extracted from the summary statistics datasets (Table 2, 3).

Two-sample MR Analysis

In the univariable MR analysis, the traditional method used was the inverse variance weighted (IVW) approach. This method combines the Wald ratio estimates for each instrumental SNP through a model similar to meta-analysis, providing an overall estimate of the effect.^[19] When there was no heterogeneity (P > 0.05), the fixed effects model (FEM) was used; if heterogeneity existed (P < 0.05), the multiplicative random-effect model (REM) was used. However, the IVW method only generates an unbiased estimate under the specific assumptions that there is no invalid instrument and horizontal pleiotropy. Several sensitivity analysis methods were employed to validate the robustness of the effect estimate obtained from the IVW method, including the weighted median (WM) and MR-Egger. The WM validly estimates the causal effect even if up to half of the weights are driven by invalid instruments.^[23] The intercept of MR-Egger regression captures the average pleiotropic effect across the genetic variants, where a significant non-zero intercept term can indicate directional pleiotropy.^[24]

Sensitivity Analysis

Cochrane’s Q value was applied to assess the heterogeneity between SNPs. Visual funnel plots were applied to evaluate publication bias and directional pleiotropy by examining symmetry.^[25] Horizontal pleiotropy was detected via conducting the MR-Egger-intercept and the MR pleiotropy residual sum and outlier (MR-PRESSO) tests. The outlier test would compare the expected and observed distribution of each variant to identify outlier variants. If any outliers were detected, they would be eliminated, and the corrected results would be re-evaluated.^[26] In addition, a leave-one-out (LOO) analysis was applied to avoid horizontal pleiotropy induced by a single SNP, which systematically eliminated the SNP and calculated the effect of the remaining SNPs.^[27] All statistical analyses were conducted using the “TwoSampleMR” package (version 0.5.7) and the “MRPRESS” package (version 1.0) in R software (version 4.3.1).

Causal Effects of OSA on Rheumatic Diseases

As mentioned above, we got 12 independent SNPs from OSA through rigorous criteria. Table S3 displays statistics for the 12 SNPs chosen as valid IVs in this MR analysis. The general and each IV’s F values were greater than the empirical threshold 10, suggesting weak instrumental factors did not cause significant bias. Then, we harmonized the summary statistics of exposure and outcome by aligning the efects to the identical reference alleles and removed abnormal, ambiguous, and palindromic SNPs from IVs. Finally, SNPs related to OSA screened in the rheumatic immune disease database were 8 SNPs in SLE and 11 SNPs in other rheumatic disease (Table 2).

There was no evidence supporting causal relationships between OSA and rheumatic disease. Using the IVW method, the MR results indicated that OSA could not directly affect the risk of rheumatic disease (all P > 0.05). These findings were confirmed using other MR methods, including WM and MR Egger (Table 2, Figure 2).

Figure 2

MR estimates for the causal effect of OSA on rheumatic disease. SNPs, single nucleotide polymorphisms; OR, odds ratio; CI, confidence interval; IVW, inverse variance weighted; FEM, fixed effects model; REM, random effects model; OSA, Obstructive sleep apnea; AS, Ankylosing spondylitis; RA, Rheumatoid arthritis; SLE, Systemic lupus erythematosus; SS, Sjogren’s syndrome; BD, Behcet disease; PsA, Psoriatic arthritis; SSc, Systemic sclerosis.

The Cochran’s Q statistic (P > 0.05), MR Egger intercept test (P > 0.05), and MRPRESSO (Global P > 0.05) failed to identify any directional pleiotropy or heterogeneity (Table S1). Besides, the forest plots, scatter plots, funnel plots, and LOO plots are displayed in Figure S1–4.

Causal Effects of Rheumatic Diseases on OSA

We conducted a reverse analysis to assess further the causal effect of rheumatic disease on OSA. After removing abnormal, ambiguous, and palindromic SNPs, we got 6 SNPs for PsA, 4 SNPs for gout, 4 SNPs for SS, 14 SNPs for BD, 8 SNPs for SSc, 48 SNPs for RA, 36 SNPs for SLE and 10 SNPs for AS. The F statistic value at each chosen IV was greater than 10, indicating that the chosen SNPs were adequately robust and that weak IVs were not likely to alter the causal estimate. The specific information can be found in the Table S6–13. Genetically predicted AS was associated with the risk of OSA (IVW: OR = 1.0239, 95% CI = 1.0086 to 1.0394, P = 0.0021; MR-Egger: OR = 1.0374, 95% CI = 1.0089 to 1.0668, P = 0.0326; weighted median: OR = 1.0287, 95% CI = 1.0109 to 1.0467, P = 0.0014). However, no evidence suggests a causal relationship between OSA and other rheumatic disease (Table 3, Figure 3).

Figure 3

MR estimates for the causal effect of rheumatic disease on OSA. SNPs, single nucleotide polymorphisms; OR, odds ratio; CI, confidence interval; IVW, inverse variance weighted; FEM, fixed effects model; REM, random effects model; OSA, Obstructive sleep apnea; AS, Ankylosing spondylitis; RA, Rheumatoid arthritis; SLE, Systemic lupus erythematosus; SS, Sjogren’s syndrome; BD, Behcet disease; PsA, Psoriatic arthritis; SSc, Systemic sclerosis.

Except for the heterogeneity of PsA, our analysis did not uncover any signs of heterogeneity or pleiotropy in the causal relationship between rheumatic disease and OSA, as assessed by Cochrane’s Q (P > 0.05) and MR Egger intercept (P > 0.05) (Table S2). Through MR-PRESSO and LOO analysis, one outlier (rs6679677) was identified in RA, and two outliers (rs6679677, rs28361029) were identified in SLE, which indicated that the IVs of exposure and outcome have significant horizontal pleiotropy. After removing outliers, the calculation was repeated. It was found that there was no horizontal pleiotropy, and these results were not significantly affected. Figure S5–8 contains displays of some plots, including forest plots, scatter plots, funnel plots, and LOO plots.