Multiple Sclerosis and Hematologic Malignancies: A Bidirectional Mendelian Randomization Study

Summary

This bidirectional two-sample Mendelian randomization study used European GWAS to assess causal links between multiple sclerosis (MS) and hematologic malignancies (HM). Genetic liability to MS was associated with higher odds of Hodgkin lymphoma and unspecified leukemia; no associations for other subtypes or in reverse analyses.

Abstract

Observational studies have reported associations between multiple sclerosis (MS) and hematologic malignancies (HM), but findings remain inconsistent. We conducted a bidirectional two-sample Mendelian randomization (MR) analysis using publicly available genome-wide association summary statistics for MS (n = 115,803) and HM (n = 218,792). Primary causal estimates were obtained with the inverse-variance-weighted (IVW) estimator, complemented by prespecified sensitivity analyses for heterogeneity and pleiotropy. Genetically proxied liability to MS was associated with higher odds of leukemia (unspecified subtype) (odds ratio [OR] 1.311, 95% confidence interval [CI] 1.002-1.716, P = 0.048) and Hodgkin lymphoma (HL) (OR 1.224, 95% CI 1.052-1.425, P = 0.009), with no evidence for other leukemia, lymphoma, or plasma-cell neoplasm subtypes (all P > 0.05). Tests indicated no substantial heterogeneity or directional pleiotropy for the leukemia (unspecified subtype) or HL analyses. These results provide genetic evidence consistent with an elevated risk of select HM subtypes among individuals with higher genetic liability to MS; however, the signals should be interpreted cautiously and validated in larger, multi-ancestry datasets and with additional causal frameworks.

Introduction

Multiple sclerosis (MS) is a chronic autoimmune disorder of the central nervous system characterized by inflammatory demyelination and a relapsing-remitting course that can lead to disability and substantial societal and familial burden¹. Hematologic malignancies (HM) comprise heterogeneous neoplasms of hematopoietic or lymphoid origin and account for approximately 6.5% of all tumors worldwide². Although targeted and immunotherapies have improved outcomes for many patients, prognosis remains poor for a subset³. The relationship between MS and HM has been investigated for decades with divergent results^4,5. In a Norwegian registry study including 6,883 MS cases and 37,919 matched controls, no statistically significant association was observed between MS and hematologic cancer risk. Similarly, a Danish study reported only a marginal-statistically non-significant-elevation in HM incidence among patients with MS compared with population controls. By contrast, some cohort studies have suggested higher risks in specific contexts or subtypes: a Korean prospective study reported a 4.52-fold higher incidence of lymphoma among individuals with MS⁶, and a large cohort analysis (29,617 MS vs. 296,164 non-MS) found higher overall HM risk in the MS group⁷. Other investigations, however, did not observe increased risks for non-Hodgkin lymphoma (NHL) or HL among patients with established MS^8,9. Collectively, the evidence for an association-particularly at the level of HM subtypes-remains uncertain.

Mendelian randomization (MR) leverages germline genetic variants as proxies for modifiable exposures under the principle that allele allocation is effectively random with respect to many environmental confounders¹⁰. By using genetic instruments strongly associated with a given exposure, MR can strengthen causal inference and reduce bias from confounding or reverse causation¹¹.

Here, we apply a bidirectional, two-sample MR framework to evaluate potential causal links between MS and the major HM classes-leukemia, lymphoma, and multiple myeloma (MM)-from a genetic perspective. Originality statement: To our knowledge, this is the first bidirectional MR to systematically interrogate subtype-specific MS-HM relationships while prespecifying operational parameters (instrument thresholds, LD clumping, harmonization, and sensitivity diagnostics) to enable full reproducibility.

Protocol

This study analyzed de-identified, summary-level genome-wide association study (GWAS) statistics that are publicly available. In accordance with repository policies and the approvals obtained by the original investigators, no new institutional review board approval or additional individual informed consent was required for this secondary analysis. All contributing GWAS reported ethics approval and consent procedures in their source publications. All analyses were conducted in compliance with institutional guidelines and the Declaration of Helsinki.

Overview and rationale

The study implemented a bidirectional, two-sample Mendelian randomization (MR) framework restricted to European-ancestry summary statistics to evaluate potential causal relationships between multiple sclerosis (MS) and hematologic malignancies (HM). The design adheres to the three core MR assumptions: instrument relevance, independence from confounders, and exclusion restriction. The workflow therefore includes (i) dataset access and curation, (ii) instrument selection at genome-wide significance with linkage disequilibrium (LD) clumping, (iii) confounder screening using PhenoScanner, (iv) allele harmonization with explicit handling of palindromic variants, (v) directionality assessment using the Steiger test¹², (vi) primary MR estimation with complementary methods, (vii) a full set of sensitivity diagnostics, and (viii) standardized figure and table generation under multiple-testing control. Each of these steps is described in detail in the subsequent protocol subsections, and an overview of the pipeline is presented in Figure 1.

Materials, software, and RRIDs

Analyses were conducted in R version 4.3.1 (RRID:SCR_001905) using RStudio/Posit 2023.12+ (RRID:SCR_000432). LD clumping, when performed locally, used PLINK v1.9 (build 2.3; RRID:SCR_001757)¹³. MR estimation and data extraction used the R package TwoSampleMR v0.5.7 ¹⁰; instrument look-ups for potential confounders used phenoscanner v1.0; detection and correction of outliers used MRPRESSO v1.0. Exact versions are reported for packages without RRIDs.

Data sources and access

MS summary statistics were obtained from the International Multiple Sclerosis Genetics Consortium meta-analysis comprising 47,429 MS cases and 68,374 controls with harmonized quality control across 15 cohorts. HM summary statistics were obtained from FinnGen (overall n = 218,792; >16 million variants) and included Hodgkin lymphoma (HL), diffuse large B-cell lymphoma (DLBCL), follicular lymphoma (FL), mature T/NK-cell lymphomas (MTNKL), other or unspecified non-Hodgkin lymphoma (NHL), lymphoid leukemia, myeloid leukemia, leukemia of unspecified cell type, and multiple myeloma/plasma-cell neoplasms¹⁴. Datasets were accessed through the IEU OpenGWAS portal using documented accession identifiers¹⁵. All analyses in this study were therefore based exclusively on these publicly available summary-level GWAS datasets; no internal institutional cohort or individual-level patient data were used or generated. Because we did not identify additional GWAS with harmonized MS and hematologic malignancy subtype definitions that would allow a full replication of the pipeline, independent external validation using a separate dataset was not performed and is acknowledged as a limitation. The protocol is written so that it can be directly re-applied to future GWAS datasets for independent validation.

Instrument selection and LD clumping

For each exposure, single-nucleotide polymorphisms (SNPs) were selected at genome-wide significance (P < 5 × 10^-8) using the extract_instruments function in TwoSampleMR applied to the OpenGWAS datasets. To ensure instrument independence, LD clumping was then performed against a European-ancestry reference panel using either the internal clumping utilities of TwoSampleMR or locally with PLINK, with an r² threshold of 0.001 and a physical window of 10,000 kilobases. When PLINK was used, the command-line parameters were set to a primary significance threshold of 5 × 10^-8, r² = 0.001, and a 10-Mb window so that the clumped instruments exactly matched these criteria. Instrument strength was evaluated using the F-statistic derived from the exposure effect estimate and its standard error (F ≈ β²/SE²); variants with F < 10 were excluded from the final instrument sets, and the remaining SNPs were carried forward to PhenoScanner screening.

Confounder screening with PhenoScanner

To minimize horizontal pleiotropy through known risk factors, each candidate instrument was queried in PhenoScanner V2 across the GWAS catalog using the phenoscanner R package (v1.0)^16,17. For every SNP, we requested all reported associations at P < 1 × 10⁻⁵ and manually inspected the returned traits. Associations indicating links to established hematologic malignancy risk factors-such as smoking-related exposure or adiposity/anthropometric traits (e.g., body mass index, waist circumference, and body fat measures)-or direct associations with hematologic malignancy phenotypes prompted exclusion of the corresponding SNP from the instrument set¹⁸. Trait categories considered as exclusionary were based on prior evidence relating obesity and smoking to leukemia, lymphoma, or myeloma risk^18,19,20. Queries used broad keyword stems (e.g., smoke, cigarette, BMI, obesity, waist, adiposity, hematologic malignancy, lymphoma, leukemia, myeloma). All removals were documented in a tracking spreadsheet along with the PhenoScanner trait that triggered exclusion, and the cleaned instrument lists were then passed to the harmonization step.

Harmonization and palindromic handling

Effect alleles for each SNP were harmonized between the exposure and outcome datasets using the harmonise_data function in the TwoSampleMR package (v0.5.7, R). We aligned all outcome alleles to the exposure effect allele so that positive beta coefficients always corresponded to the same allele in both datasets. Palindromic variants (A/T or C/G) with intermediate effect-allele frequencies (0.42-0.58) in the OpenGWAS reference panel were treated as strand-ambiguous and automatically removed by setting the harmonization action to drop ambiguous SNPs. Palindromic SNPs with effect-allele frequencies outside this range were retained and aligned using the reported allele frequencies. Because allele availability and palindromic status differed slightly across FinnGen outcomes, harmonization was run separately for each HM phenotype, and the final number of instruments entering each outcome-specific analysis was extracted from the harmonized R objects and reported in the tables.

Directionality assessment (Steiger filtering)

Directionality was evaluated using the Steiger approach as implemented in the steiger_filtering function of TwoSampleMR. For each SNP, the function first calculated the variance explained (R²) in the exposure and outcome from the GWAS beta coefficient, standard error, and sample size. The study then removed instruments for which R² was greater in the outcome than in the exposure, indicating a possible reverse direction of effect. Steiger filtering was applied separately for each outcome dataset, and the remaining instruments (rows with steiger_dir == TRUE) were saved and used in the subsequent MR analyses. Post-Steiger instrument counts were recorded for each outcome and are reported alongside the MR estimates.

Primary MR estimation and multiple-testing control

Primary causal estimates were obtained with inverse-variance-weighted (IVW) MR under a fixed-effects model using the mr function in TwoSampleMR, with methods specified as "mr_ivw", "mr_egger_regression", and "mr_weighted_median". For each HM outcome, harmonized and Steiger-filtered instruments were passed to mr, and log odds ratios and standard errors were extracted and exponentiated to obtain odds ratios (ORs) with 95% confidence intervals (CIs) for binary traits²¹. To examine robustness to modest violations of the no-pleiotropy assumption, we additionally applied the Weighted Median and MR-Egger regression estimators^22,23, implemented in the same package. When Cochran's Q test (from mr_heterogeneity) indicated substantial heterogeneity (P < 0.05), the study also fitted multiplicative random-effects IVW models and reported both fixed- and random-effects results. Family-wise error across the nine HM outcomes was controlled using Bonferroni correction with α = 0.05/9 = 5.56 × 10^-3; associations with P values below this threshold were considered statistically significant, whereas those with 0.0056 ≤ P < 0.05 were interpreted as suggestive and described cautiously.

Sensitivity diagnostics: Heterogeneity, pleiotropy, and outliers

Cochran's Q statistic was used to assess between-instrument heterogeneity for both IVW and MR-Egger models, implemented via the mr_heterogeneity function in TwoSampleMR. Directional horizontal pleiotropy was evaluated using the MR-Egger intercept test (mr_pleiotropy_test) and the global test in the MR-PRESSO package²⁴. MR-PRESSO²⁴ was run with the recommended settings in R (NbDistribution ≥ 5,000, SignifThreshold = 0.05) to detect influential outliers and to quantify potential distortion by comparing IVW estimates before and after outlier removal²⁵. Leave-one-out analyses (mr_leaveoneout) were performed for each exposure-outcome pair to determine whether any single SNP disproportionately influenced the overall estimate. For transparency and reproducibility, all diagnostic outputs were exported from R and reported together with the corresponding instrument counts after harmonization, Steiger filtering, and MR-PRESSO outlier removal.

Instrument strength and NOME assessment

Instrument strength for MR-Egger was quantified using the I²GX statistic, computed as 1 minus the mean of the squared standard errors of the SNP-exposure associations divided by their variance across instruments²⁶. Values closer to 1 indicate better compliance with the No Measurement Error (NOME) assumption; lower values suggest possible regression dilution and prompt cautious interpretation of MR-Egger results. I²GX was calculated and reported for each outcome-specific analysis.

Reverse Mendelian randomization

The complete pipeline was repeated in the reverse direction by treating each HM subtype as the exposure and MS as the outcome. When genome-wide significant instruments were insufficient for a given HM exposure, a relaxed selection threshold of P < 5 × 10^-6 was allowed while maintaining the same LD clumping parameters, PhenoScanner screening, harmonization procedures, Steiger filtering, and sensitivity diagnostics. Analyses that used relaxed thresholds were clearly labeled in the corresponding tables and figure legends.

Visualization and figure export

Scatter, forest, funnel, and leave-one-out plots were generated with legends positioned below the panels and font sizes adjusted to ensure that labels did not obscure plotted data. Axis limits were standardized across comparable outcomes to facilitate visual comparison. Figures were exported at a minimum of 300 dpi in lossless formats such as TIFF or PNG. All plotted numerical values were cross-checked against the reported estimates to ensure consistency between text, tables, and figures.

Reproducibility and data sharing

Random seeds were fixed where applicable, software versions were recorded, and analysis scripts together with intermediate objects were archived to allow re-running of all steps. Dataset accession identifiers and phenotype definitions were documented, and the instrument lists at each filtering stage-post-clumping, post-harmonization, post-Steiger filtering, and post-MR-PRESSO were prepared for upload as spreadsheet files in accordance with journal guidelines.

Results

Selection of instrumental variants

Using a genome-wide significance threshold of P < 5×10⁻⁸, 22,466 SNPs were initially retrieved for MS. After LD clumping in PLINK (r² < 0.001, window 10,000 kb), 72 SNPs remained as instruments, and each had F > 10 (range 30-1,044). To minimize pleiotropy through known risk pathways, 11 SNPs associated with potential HM confounders (e.g., body mass index, waist circumference, weight, hip circumference, fat percentage, fat ...

Discussion

This two-sample Mendelian randomization (MR) study in a European-ancestry population provides genetic evidence consistent with a higher risk of Hodgkin lymphoma (HL) and leukemia of unspecified cell type among individuals with greater genetic liability to multiple sclerosis (MS), while showing no convincing association with other leukemia subtypes, non-Hodgkin lymphoma (NHL) subtypes, or multiple myeloma/plasma-cell neoplasms²⁷. Signals for HL and unspecified leukemia are suggestive rather than de...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
Bonferroni Correction (α = 0.0056)	Statistical method reference	BC-0056	Multiple testing correction method to control the false positive rate
Cochran's Q Statistic	Cochran, 1954	Q-001 (https://www.jstor.org/stable/2334368)	Used to assess heterogeneity between MR instruments
F-statistic	TwoSampleMR documentation	F-049 (https://cran.r-project.org/web/packages/TwoSampleMR/vignettes/perform_mr.html)	Metric for assessing instrument strength in MR
HM Summary Statistics (GWAS)	FinnGen	https://www.finngen.fi/en	Genomic data for hematologic malignancies (n = 218,792), including Hodgkin lymphoma (HL), non-Hodgkin lymphoma (NHL), leukemia, etc.
I2GX statistic	Bowden et al., 2016	I2GX-074 (https://academic.oup.com/ije/article/44/2/512/753845)	Statistic to evaluate compliance with the "no measurement error" assumption
LD Clumping Parameters	PLINK	https://www.cog-genomics.org/plink/1.9	Parameters for LD clumping to ensure independence of genetic instruments
MR-PRESSO v1.0	MR-PRESSO	1445 (https://github.com/rondolab/MR-PRESSO)	Software for heterogeneity and outlier detection and correction
MS Summary Statistics (GWAS)	International MS Genetics Consortium	https://imsgc.net	Genomic data for multiple sclerosis (MS) (n = 115,803), from European cohort
PhenoScanner V2	PhenoScanner	1550 (http://www.phenoscanner.medschl.cam.ac.uk)	Database to screen genetic instruments for associations with known confounders (e.g., smoking, BMI)
PLINK v1.9	PLINK	SCR_001757 (https://www.cog-genomics.org/plink/1.9)	Software for genome-wide association analysis and linkage disequilibrium (LD) clumping
R version 4.3.1	R Foundation	SCR_001905 ( https://www.r-project.org)	Open-source statistical computing software
Steiger Test	TwoSampleMR methodology	STG-001 (https://cran.r-project.org/web/packages/TwoSampleMR)	Statistical test for directionality in Mendelian randomization
TwoSampleMR v0.5.7	RStudio / CRAN	1134 (https://cran.r-project.org/package=TwoSampleMR)	R package for two-sample Mendelian randomization (MR) analysis
Weighted Median and MR-Egger	MendelianRandomization methods	WM-026 ( https://cran.r-project.org/package=MendelianRandomization)	Alternative estimators for sensitivity analysis in MR