探究腰痛、髋关节痛、膝关节痛与步态异常之间的双向关系：一项孟德尔随机化研究

摘要

IVW 分析表明，下背部和髋部疼痛可能与步态异常有关，而膝关节疼痛及反向关联尚无定论。

摘要

步态异常与下肢生物力学改变及肌肉骨骼系统疾病中的功能障碍相关。尽管观察性研究已将腰痛、髋部疼痛和膝关节疼痛与步态受损联系起来，但它们之间的因果关系仍不明确。我们开展了一项双向孟德尔随机化（MR）研究，以评估特定部位疼痛与步态异常之间潜在的因果关联。采用来自大规模全基因组关联研究的独立单核苷酸多态性（SNPs）作为腰痛、髋部疼痛和膝关节疼痛的遗传工具变量。在正向MR分析中，将源GWAS中定义为自我报告的行走困难的步态异常作为结局指标；在反向MR分析中则作为暴露因素。主要分析采用逆方差加权法（IVW），并辅以敏感性分析以评估异质性、水平多效性和结果稳健性。在正向分析中，基于遗传预测的腰痛（OR = 1.53，95% CI = 1.058–2.219，P = 0.024）和髋部疼痛（OR = 1.55，95% CI = 1.067–2.252，P = 0.021）在IVW分析中与步态异常风险增加相关。然而，其余四种MR方法的结果均无统计学显著性，尽管效应方向总体一致。遗传预测的膝关节疼痛未显示显著关联（OR = 2.15，95% CI = 0.234–19.789，P = 0.498），其宽置信区间反映了较大的不确定性。反向MR分析未发现步态异常的遗传易感性会增加腰痛、髋部疼痛或膝关节疼痛风险的证据，但该分析仅能利用四个与步态异常相关的SNP。未检测到显著的异质性或水平多效性。总体而言，本研究提供了有限的基于IVW的遗传学证据，提示腰痛和髋部疼痛可能与步态异常相关。研究结果应谨慎解读，并需利用更大样本量的GWAS数据集及更精确的表型进行验证。

引言

步态异常可能源于运动或感觉功能障碍，其临床特征取决于潜在病理的部位和性质。特定的异常步态模式可为某些疾病提供重要的诊断线索。值得注意的是，在肌肉骨骼系统疾病中，步态异常常表现出独特的生物力学特征，并与疼痛、关节功能障碍及活动能力下降密切相关^1,2,3。研究表明，长期的步态异常可加重骨关节炎患者的病情⁴。此外，持续的异常步态可能增加运动关节的不稳定性并改变关节负荷；在异常生物力学条件下反复进行的关节活动可能进一步加剧关节损伤，从而引发疼痛，促进肌肉逐渐萎缩，并导致更明显的肌力丧失。这些因素可能相互作用，形成恶性循环，促使肌肉骨骼系统疾病患者的疼痛和功能障碍进行性加重⁵。因此，明确特定部位肌肉骨骼疼痛与步态异常之间关联的方向性及潜在因果关系，有助于解释观察到的疼痛–步态关系，并为未来关于步态功能障碍及康复的研究提供依据。

然而，尽管越来越多的观察性证据表明特定部位的疼痛（如腰痛、髋关节疼痛和膝关节疼痛）与步态功能受损相关，但传统的观察性研究容易受到混杂因素、测量误差和反向因果关系的影响，使得难以确定这些关联的方向和因果关系。此外，疼痛与步态之间的关系可能涉及中枢神经系统的适应以及来自外周的伤害性信息的调节，这进一步增加了对所观察到的疼痛-步态关联解释的复杂性。孟德尔随机化（Mendelian randomization, MR）是一种旨在克服观察性研究关键局限性的流行病学方法，已广泛应用于因果推断研究中。MR利用独立的单核苷酸多态性（SNPs）作为工具变量，推断暴露因素与结局之间的潜在因果关系⁶。例如，MR无需人为让个体产生腰痛，而是可以检验在大规模全基因组关联研究（GWAS）数据集中，个体在遗传上易患腰痛的倾向是否与步态相关结局存在关联。由于基因变异在受精时即被随机分配，即在胚胎形成时确定，因此MR可减少年龄、体重指数、体力活动以及共病性肌肉骨骼疾病等混杂因素带来的偏倚。当具备统计效能充足的GWAS数据和有效的遗传工具变量时，该方法适用于探索人群水平的因果关系。然而，本研究中使用的表型为广泛的自我报告型GWAS性状，可能无法充分反映疼痛严重程度、症状持续时间、经临床医生确认的诊断或客观的步态参数。通过利用基因变异在受精时的随机分配特性，MR能够有效减轻混杂偏倚和反向因果关系所带来的影响，从而提供更为可信的因果证据。

本研究的目的是利用关于自报腰痛、髋关节痛和膝关节痛的全基因组关联研究（GWAS）汇总数据，结合IEU Open GWAS数据库中关于步态异常的GWAS数据，进行双向孟德尔随机化（MR）分析。在本研究中，疼痛表型代表了广泛的自报部位特异性疼痛特征，而步态异常则指自报行走困难，而非临床医生确诊的诊断、实验室测量的步态参数或标准化的功能评分。该双向设计旨在利用GWAS汇总统计量评估部位特异性肌肉骨骼疼痛与步态异常之间潜在因果关系的方向⁷。

方案

本研究使用了来自 IEU OpenGWAS 数据库和 GWAS 目录的公开可用的全基因组关联研究（GWAS）汇总级别数据。未访问任何个体层面的参与者数据，也未招募新的研究对象。原始贡献性 GWAS 研究中已获得伦理批准和书面知情同意。因此，本研究对公开可用的汇总统计资料进行二次分析无需额外的机构审查委员会批准。

研究原理与设计

MR 研究的基本步骤包括获取 GWAS 汇总数据、筛选和评估 SNPs、进行统计分析以及实施质量控制措施。如图1所示，MR 分析的准确性依赖于满足三个关键假设：（1）相关性假设：工具变量（IVs）必须与暴露表型相关，即腰痛、髋关节痛和膝关节痛⁸；（2）独立性假设：工具变量（IVs）与影响"暴露-结局"的混杂因素无关；（3）排他性假设：工具变量（IVs）仅通过暴露途径影响结局，而不通过其他路径⁹。在每次分析中，暴露因素和结局变量将互换位置，以判断两者之间是否存在反向因果关系。

为了提高可重复性，用于本地GWAS汇总统计量导入、SNP提取、LD剪枝、结局SNP提取、等位基因匹配、MR分析、敏感性分析、诊断绘图以及结果保存的完整计算工作流程均作为补充代码1提供。正文部分描述了关键的实验步骤，而补充文件1则提供了相应的R语言命令级和函数级实现方法。

数据来源

与腰痛、髋关节疼痛和膝关节疼痛相关的数据来自 IEU openGWAS 数据库，网站地址为：https://gwas.mrcieu.ac.uk/。用于定义这些疾病的具体问题见补充文件 2。在原始 GWAS 数据集中，腰痛、髋关节疼痛和膝关节疼痛是根据自我报告的部位特异性疼痛项目来定义的。腰痛数据集的 GWAS ID 编号为：ebi-a-GCST90018797，样本量为 468269，其中包括 22,413 名腰痛患者和 445,856 名对照个体，总共包含 24,174,741 个 SNP。髋关节疼痛数据集的 GWAS ID 编号为：ebi-a-GCST90013968，样本量为 407,746，SNP 总数为 11,039,206。膝关节疼痛数据集的 GWAS ID 编号为：ukb-b-16254，样本量为 461,857，其中包括 98,704 名膝关节疼痛患者和 363,153 名对照个体，总共包含 9,851,867 个 SNP。步态异常的 GWAS 数据也来自 IEU openGWAS 数据库，步态异常数据集的 GWAS ID 编号为：finn-b-R18 _ABNORMALITI_GAIT_MOBIL，样本量为 210,717，包括 1,348 例步态异常病例和 209,369 名对照个体，共鉴定出 1,638,044 个 SNP 位点。详细信息见表 1。

为确保数据获取的可重复性，研究人员应在 IEU OpenGWAS 数据库中逐一查询上述每个 GWAS ID，下载相应的 GWAS 汇总统计文件（如可获取），并将文件保存为本地文本文件或逗号分隔文件，然后再导入 R 中。在本研究采用的计算流程中，根据源文件格式，使用 read.table() 或 read.csv() 函数将本地 GWAS 汇总统计文件导入 R。采用基于本地文件的工作流程，是为了确保相同的下载 GWAS 汇总统计文件能够通过固定的列映射和完全一致的分析命令进行重复处理。导入和处理本地文件所使用的具体 R 命令详见补充文件 1。

SNP的筛选与录入

使用全基因组显著性阈值 P < 5 × 10⁻⁸，从 GWAS 汇总数据集中提取与每种暴露表型显著相关的单核苷酸多态性（SNP）。若在此阈值下可用的独立 SNP 少于三个，则将阈值放宽至 P < 5 × 10⁻⁶，以获得足够的遗传工具用于孟德尔随机化（MR）分析。通过在 10,000 kb 窗口内使用 r² < 0.001 评估连锁不平衡，以确保工具变量的独立性。

在脚本层面，首先将暴露性状的GWAS汇总统计结果导入R，并根据预设的P值阈值进行筛选。筛选出的暴露性状相关SNP被保存为一个暴露文件，并按照TwoSampleMR软件包的要求进行格式化。所需的列映射包括SNP标识符、效应估计值、标准误、效应等位基因、其他等位基因以及P值。在可重复的R工作流程中，使用read_exposure_data()函数导入暴露数据，并通过设置clump = TRUE执行连锁不平衡（LD）剪枝。相应的R命令结构为：read_exposure_data(filename = "exposure.csv", sep = ",", snp_col = "rsids", beta_col = "beta", se_col = "sebeta", effect_allele_col = "alt", other_allele_col = "ref", pval_col = "pval", clump = TRUE)。列名称根据实际下载的GWAS文件进行相应调整。

使用 F 统计量评估工具变量的强度：

F = [(N − k − 1)/k] × [R²/(1 − R²)]，

其中，N 表示样本量，k 表示工具变量的个数，R² 表示SNP所解释的暴露方差比例。R² 被计算为

Σ[2 × MAF × (1 − MAF) × β²/(SE² × N)]，

其中，MAF 为次要等位基因频率，β 为等位基因效应估计值，SE 为标准误。

保留F统计量>10的SNP。使用PhenoScanner V2识别并剔除与潜在混杂因素（包括先天性解剖异常和体重指数）相关的SNP。

通过将 clumped 的暴露 SNP 与相应的局部结局 GWAS 汇总统计文件根据 SNP 标识符进行合并，提取结局 SNP。如果结局 GWAS 报告的是 −log10(P)，则使用 P = 10^(−LP) 转换为 P 值。提取的结局 SNP 随后通过 read_outcome_data() 导入。使用 harmonise_data() 对暴露和结局数据集进行整合，以统一效应等位基因并剔除不适用于 MR 分析的 SNP。具有模糊等位基因方向的回文 SNP 被排除，而 mr_keep = TRUE 的 SNP 被保留用于后续的 MR 分析¹¹。结局 SNP 提取、P 值转换和数据整合的详细脚本见补充文件 1。

统计分析

本研究采用五种方法估计暴露变量与结局变量之间的因果效应：逆方差加权法（IVW）、MR-Egger回归法、加权中位数法、加权众数法和简单众数法¹²。IVW法被认为是孟德尔随机化（MR）分析的标准方法，其假设所有工具变量均为有效变量，基于合并Wald比率估计值的原理；若存在异质性，则采用随机效应模型，若无异质性，则采用固定效应模型。MR-Egger回归法可检测潜在的多重共线性，并在回归中考虑截距项的存在¹³。加权中位数法要求超过50%的工具变量为有效的SNP。加权众数法相比其他方法所需样本量更小，且能减少偏倚并降低Ⅰ类错误率。简单众数法允许根据所估计的因果效应是否相似，将具有相似效应的SNP进行分组¹⁴。

在功能层面，使用 TwoSampleMR 软件包中的 mr() 函数进行 MR 分析。分析中使用的方法列表指定为 c("mr_ivw", "mr_egger_regression", "mr_weighted_median", "mr_simple_mode", "mr_weighted_mode")。比值比及其 95% 置信区间通过 generate_odds_ratios() 函数生成。完整的 R 语言实现代码，包括确切的方法列表和输出保存命令，详见补充文件 1。

质量控制与敏感性分析

为了检验 MR 结果的稳定性和可靠性，首先采用 Cochran Q 检验评估 SNP 之间的异质性，若 Cochran Q 检验具有统计学显著性，则表明分析结果中存在显著的异质性¹⁵。其次，采用 MR-Egger 截距检验评估潜在的水平多效性。若 MR-Egger 截距具有统计学显著性，则提示 MR 分析中存在方向性水平多效性。第三，应用孟德尔随机化多效性残差和离群值（MR-PRESSO）方法检测结果中是否存在离群 SNP，若存在，则将其剔除后重新分析¹⁶。第四，采用"逐一剔除"方法检验结果的稳健性，通过逐个排除 SNP 并计算剩余 SNP 的合并效应，以评估单个 SNP 对暴露与结局变量之间关联的影响¹⁷。

在功能层面，使用 mr_heterogeneity() 评估异质性，使用 mr_pleiotropy_test() 评估方向性多效性，使用 mr_singlesnp() 评估单个 SNP 效应，并使用 mr_leaveoneout() 进行留一法分析。诊断图通过 mr_scatter_plot()、mr_funnel_plot() 和 mr_leaveoneout_plot() 生成。当可用工具变量数量足够时，执行 MR-PRESSO 分析；若因 SNP 数量不足而无法执行 MR-PRESSO，则记录并报告该情况。MR 分析及质量控制过程使用 R 4.3.2 版本和 TwoSampleMR 软件包 0.5.8 版本完成，显著性水平设定为 α = 0.05。完整命令见补充文件 1。

输出保存、质量检查与报告

所有中间和最终输出结果均被保存，以确保可重复性。这些结果包括：选定的暴露相关SNP、经过clumping处理的暴露工具变量、提取的结局相关SNP、经过数据协调（harmonized）的数据集、保留的SNP列表、孟德尔随机化（MR）估计值、比值比、异质性检验结果、MR-Egger截距项结果、MR-PRESSO结果、单个SNP分析结果、留一法分析结果、诊断图以及R会话信息。在报告前，对保留的SNP进行了以下检查：全基因组显著性、LD独立性、等位基因协调状态、模糊的回文变异位点、工具变量强度、异质性、定向水平多效性以及离群值SNP。主要的MR结果以比值比及其95%置信区间和P值形式报告。敏感性分析结果和诊断图在补充表格和图中提供。保存输出结果和检查文件的命令详见补充文件1。

结果

本研究分析了三种部位特异性肌肉骨骼疼痛表型（包括腰痛、髋关节疼痛和膝关节疼痛）以及一种步态异常表型的全基因组关联研究（GWAS）数据。每个暴露因素的工具变量（IVs）详细信息见补充文件3。

基于磁共振分析特定部位肌肉骨骼疼痛对步态异常的影响

本研究排除了具有连锁不平衡和回文结构的SNP，并剔除了与混杂因素相关的SNP，最终纳入的SNP将作为孟德尔随机化（MR）分析的工具变量。所有SNP的F统计量均大于10（20.86–91.68），表明弱工具变量偏倚的可能性较低。在主要的IVW分析中，基因预测的腰痛和髋部疼痛与步态异常风险升高相关。具体而言，IVW估计值在腰痛（OR = 1.53，95% CI = 1.058–2.219，P = 0.024）和髋部疼痛（OR = 1.55，95% CI = 1.067–2.252，P = 0.021）中均具有统计学意义。在主要的IVW分析中，基因预测的腰痛和髋部疼痛与步态异常风险升高相关。具...

讨论

步态异常可改变下肢关节的负荷和运动策略，可能加速骨关节炎进程，并导致功能受限，显著影响生活质量。这一临床相关性使得步态功能障碍与肌肉骨骼疼痛之间的相互作用受到越来越多的关注。观察性研究表明，特定部位的肌肉骨骼疼痛可能与步态受损有关；例如，一项针对年龄在65岁及以上的社区老年人的横断面研究发现，多部位疼痛与较慢的步速相关（n = 176）¹⁸。然而，疼痛与步态之间关系的方向性和因果性仍难以确定，因为随机对照试验（RCTs）往往难以实施，而观察性研究结果易受混杂因素和反向因果关系的影响，这可能阻碍基于证据的疼痛管理和康复治疗。

近期多项研究已利用孟德尔随机化（MR）方法探索特定部位肌肉骨骼疼痛的风险或保护因素^19,20,21。当随机对照试验（RCT）不可行或不符合伦理时，MR通过使用遗传变异作为工具变量，为加强因果推断提供了一种有效的方法

材料

本文使用的材料清单

姓名	公司	目录编号	评论
IEU OpenGWAS 数据库	MRC 综合流行病学单位，布里斯托大学	N/A	暴露和结局表型的全基因组关联研究（GWAS）汇总统计量来源。
MR-PRESSO R 软件包	Marie Verbanck / R 软件包	N/A	用于检测异常值并评估水平多效性。
PhenoScanner V2	剑桥大学	V2	网络资源，用于识别与潜在混杂因素相关的单核苷酸多态性（SNP）。
R	统计计算 R 基金会	版本 4.3.2	用于孟德尔随机化分析和质量控制的统计计算环境。
TwoSampleMR R 软件包	MRC 综合流行病学单位，布里斯托大学	版本 0.5.8	用于进行逆方差加权法（IVW）、MR-Egger、加权中位数、加权众数和简单众数的孟德尔随机化分析。