基于孟德尔随机化研究类风湿性关节炎与骨质疏松症的因果关系

朱凯 刘宛欣 王秋根

上海中医药大学附属岳阳中西医结合医院,上海 200437

骨质疏松症(osteoporosis,OP)是一种系统性骨骼疾病,以骨量丢失和骨组织微结构退化为主要特征,从而降低骨骼的机械强度,增加骨折的风险[1]。有研究发现全球OP患病率约为18.3%[2],而我国OP患病率可达20%,在女性群体中患病率更高[3]。此外,骨质疏松性骨折是老年人发病和残疾的主要原因之一,在髋部骨折的情况下甚至可导致死亡。因此,骨质疏松症正逐渐成为人口老龄化的主要健康问题[4]。类风湿性关节炎(rheumatoid arthritis,RA)是一种慢性炎症性自身免疫性疾病,其主要特征是持续性的滑膜炎及多关节骨骼和软骨的进行性破坏。流行病学研究显示全球RA的患病率为0.3%～1%[5],然而一项基于全球57项研究的Meta分析表明,OP在全球RA患者中的患病率为27.6%,其中亚洲的患病率为30.6%[6]。另一项基于44项研究的Meta分析显示,在中国RA患者中,OP的总体患病率为37.67%,远高于亚洲及全球平均水平。目前许多观察性研究表明OP与RA关系密切[7-8],然而由于反向因果关系和潜在的混杂因素影响,其因果推论容易产生偏差。此外,随机对照试验是流行病学研究中因果推论的金标准,但由于医学伦理、研究时间及成本等原因,一些大型及高质量随机对照试验也难以开展。

孟德尔随机化(mendelian randomization,MR)是一种较为可靠的研究方法,由Katan于1986年首次提出。其基于孟德尔第二遗传定律,利用遗传变异作为工具变量确定暴露因素和结局之间的因果关系,可以有效地消除混杂因素的干扰[9]。而双样本MR分析通过使用全基因组关联研究(genome-wide association studies,GWAS)公布的汇总数据大大增加了MR的范围和统计能力[10-11]。因此,本研究进行了两组样本汇总数据MR分析,以评估RA与OP的因果关系。

1 资料与方法

1.1 研究设计

本研究将RA作为暴露因素,与RA显著相关的单核苷酸多态性(single nucleotide polymorphism,SNP)作为工具变量,OP则作为结局变量。使用两样本 MR 分析的方法进行因果关联分析,利用 Cochran Q 检验进行异质性检验,最后对所得结果进行敏感性分析以验证其可靠性。为了筛选出可靠的工具变量进行两样本MR分析,本文采用了以下3个关键假设:①工具变量与RA之间存在显著的关联;②工具变量与RA-OP关联的所有混杂因素均不相关;③工具变量只能通过与RA的关联影响结局。

1.2 数据来源

本研究从(https://gwas.mrcieu.ac.uk)公共网站获取关于RA和OP的GWAS数据。RA的GWAS数据来自欧洲生物信息研究中心,ID为ieu-a-832,由Okada等[12]在2014年发表。Zhang等[13]使用该GWAS数据探究了RA与血脂异常的因果关系,此外,Yuan等[14]也利用此数据探讨了RA与冠心病及卒中的因果关联。其中共涉及58 284名欧洲男女后裔(包括14 361例RA患者和43 923例对照),SNPs 数量为8 747 963。OP的GWAS数据则来自于芬兰数据库(The FinnGen Biobank),ID为finn-b-M13_OSTEOPOROSIS,其中包括 212 778名欧洲男女后裔(包括3 203例OP患者和209 575例对照),SNPs 数量为16 380 452。研究人群的遗传背景均来自于欧洲人,以排除因种族相关混杂因素而造成的偏差。

1.3 工具变量

本研究设置全基因组显著性水平P<5.0×10-8,同时为排除连锁不平衡(linkage disequilibrium,LD)偏倚问题,将与RA相关的显著SNPs的LD参数设置如下:R2<0.001、遗传距离为10000 kb,得到与RA显著相关的有效SNPs作为工具变量,记录效应等位基因(effect allele,EA)、其他等位基因(other allele,OA)、等位基因效应大小(β)、标准误(Se)和P值等相关信息。使用PhenoScanner v2数据库(http://www.phenoscanner.medschl.cam.ac.uk/)对混杂SNPs(如与BMI、炎症性肠病、糖尿病相关)逐个查找并剔除,以满足第二假设的要求[15]。此外,SNPs对暴露的影响应对应于同一等位基因对结局的影响,故删除工具变量中的回文SNPs。为了评估IVs的强度,本研究使用以下公式计算每个SNP的F统计量。F=R2(N-K-1)/[K(1-R2)],其中N是样本量,K是工具变量的数量,R2是每个工具变量所解释的暴露变异程度,由公式R2=2×EAF×(1-EAF)×(β/SD)2计算得到[16]。如果F>10,则认为IVs与暴露之间的相关性较强,MR分析结果可以避免受到弱工具偏差的影响[17]。

1.4 孟德尔随机化分析

本研究使用了 R version 4.2.1 中的“TwoSampleMR”包进行 MR分析。采用3种不同的孟德尔随机化分析方法来评估RA对OP的影响:逆方差加权法(inverse variance weighted,IVW)、MR-Egger 回归法(MR-Egger regression,MER)、加权中位数法(weighted median,WME)。IVW使用Wald比值方法计算出每个SNPs对应的暴露结果效应,然后进行加权线性回归。当IVs满足3个假设条件时,IVW方法通常会具有较高的准确性[18]。因此,IVW常作为首选方法。为避免未知和不可测量的混杂因素的干扰,本研究将采用 MER、WME分析方法对得出的 MR结果进行检验以辅助分析。MER方法可以不受IVs有效性的影响而提供相对稳健的估计,并通过回归斜率和截距给出现有水平多效性的调整结果[19-20],但这种方法产生的评估精度相对较低[20]。若假设至少50%的IVs都是有效的,则加权中位数方法可以提供准确和稳健的效果评估[21]。最后,使用在线工具(https://shiny.cnsgenomics.com/mRnd/)进行统计学效能分析[22]。

1.5 敏感性分析

使用MR-Egger回归法和漏斗图检测多效性。MR-Egger方法可以评估和调整所选遗传变异的水平多效性的影响[19]。如果漏斗图存在不对称性,也可以检验水平多效性。其次,应用Cochrane Q统计法以及留一法来评估SNPs之间的异质性。在Cochrane Q统计中,若P>0.05,则表示没有异质性。此外,在留一法分析中,通过依次剔除每个SNPs,然后对其余的SNPs进行IVW分析,以验证是否存在对结果有强烈影响的异常值。最后,提供散点图和森林图来显示MR分析的结果。

2 结果

2.1 工具变量

本研究在排除连锁不平衡后,获得46个SNPs,剔除rs6679677、rs62395855、rs12539741、rs10790268、rs212389共5个混杂SNPs。此外,经查找存在以下6个回文SNPs:rs13330176、rs225433、rs2661798、rs34536443、rs3799963、rs4452313。删除后共得到35个SNPs,作为本研究的工作变量,详见表1。在本研究中,单个SNPs对应的F值分布范围为29.84～1143.80,F值均大于10,表明结果不存在弱IVs偏倚,说明本研究结果可靠。

表1 工具变量具体信息Table 1 Details of instrumental variables

2.2 孟德尔随机化分析结果

孟德尔随机化结果支持RA遗传易感性与OP风险增加之间的因果关系,且基于58 284的样本量和5%的alpha值,MR分析有80%的统计学效能。IVW法、MER法及WME法分析结果方向一致。其中IVW结果(OR=1.11,95%CI:1.04～1.18,P<0.001)表明,类风湿性关节炎与骨质疏松症的患病风险呈正相关,换言之,将类风湿性关节炎的基因数据对数变换后,每降低一个标准差,则会导致骨质疏松症的患病风险增加10%,且差异具有统计学意义。此外,WME结果也进一步证明类风湿性关节炎会增加骨质疏松症的患病风险(OR=1.10,95%CI:1.03～1.18,P=0.007),MR分析结果及森林图见图1～3。

图1 三种不同MR方法的分析结果Fig.1 Analysis results from the three different MR methods

图2 RA对OP因果影响的散点图Fig.2 Scatter plot about the causal effect of RA on OP

图3 RA对OP因果影响的森林图Fig.3 Forest plot about the causal effect of RA on OP

2.3 敏感性分析结果

经 Cochran Q 检验,IVW(Q=45.801,P=0.085)及 MR-Egger(Q=44.830,P=0.082)的P值均大于0.05,表明不存在异质性。此外,在多效性检验中,MR-Egger回归结果为(egger_intercept=0.008,se=0.009,P>0.05),表明本研究结果不存在水平多效性。漏斗图结果显示当逐个以 SNPs 为 IV 时,因果关联效应的散点基本呈对称分布(图4),提示不存在水平多效性,进一步验证了MR-Egger回归结果。留一法检验结果显示依次剔除各SNPs后,剩余SNPs的IVW分析结果与纳入全部SNPs的分析结果相近,未发现对RA和OP因果关联估计值影响较大的异常SNP(图5),提示结果不存在异质性,具有稳健性。

图4 漏斗图Fig.4 MR funnel plot

图5 “留一法”分析图Fig.5 Plot of leave-one-out analysis

3 讨论

本研究基于开放数据库的遗传数据,对RA是否与罹患OP风险相关进行了MR分析。为避免来自更多异质人群的潜在混淆,笔者将MR研究人群的遗传背景限制在欧洲血统的群体。MR分析显示,IVW结果与WME结果均表明RA与OP风险呈正相关,即RA患者的OP患病风险将会增加。此外,本研究中多种敏感性分析也保证了结果的稳定性和准确性。

本研究通过MR研究发现,RA患者的遗传易感性大约会增加10%OP的患病风险,这与以往一些临床研究的结果相一致。例如,越南一项队列研究显示早期RA患者所有部位骨质疏松的发生率显著高于对照组[23]。同样,在英国一项纳入304例RA患者的队列研究中,OP的总体患病率为26.5%,明显高于由性别和年龄匹配的对照组[24]。此外,日本一项临床观察研究也发现RA与OP具有相关性[25],且继发于RA的OP患者的脆性骨折发生率约为普通人群的2倍[26]。韩国一项纳入299例RA患者的回顾性研究也表明,其OP的患病率约为健康受试者的1.9倍[27]。然而,Liu等[28]研究认为没有证据表明RA与OP之间存在因果关联,这与本研究结果不同,可能是由于他们研究的结局来自多部位骨密度的GWAS数据,而不是OP的GWAS数据。因此,基于上述研究结果,笔者仍认为RA与OP之间存在因果关联。

目前关于RA增加OP患病风险的机制存在多种学说,主要涉及炎症状态、细胞因子释放和自身抗体等。慢性炎症是RA的一个重要特征,有研究发现,RA的炎症反应主要是由TNF-α、IL-6、IL-1等细胞因子分泌增多驱动的,这些细胞因子可以直接或间接激活破骨细胞,导致骨丢失;同时炎症细胞因子亦可以阻止成骨细胞的分化[29]。也有研究表明,IL-34可促进浸润关节的炎性单核细胞分化为破骨细胞,加快骨吸收[30]。此外,炎症可通过全身或局部释放能够直接降解骨组织的蛋白酶(如金属蛋白酶),进而导致骨质疏松症[29]。另有研究发现,RANKL是参与RA局部和全身性骨丢失发病机制的关键细胞因子之一,在RA患者中,RANKL的主要来源于CD4+及CD28-T细胞,在这种情况下,RANKL可增强破骨细胞活性,同时对成骨细胞的发育产生抑制作用[31-32]。Dkk-1是Wnt信号抑制剂,是骨重塑的另一个主要调节因子,其水平增加与RA患者罹患OP和骨侵蚀的风险相关[33-34]。也有证据支持抗瓜氨酸蛋白抗体可直接激活破骨细胞[35]。尽管目前机制十分复杂,但充分说明了两种疾病在病理生理学等方面具有共同特性,最终这些机制都会影响到成骨细胞和破骨细胞的动态平衡,进而影响骨量导致OP的发生。

本研究采用两样本MR研究方法,从遗传学角度评估了RA对OP的因果效应,具有一定的优势。首先,MR分析中的遗传变异在受孕时就已分配,不受外部因素的干扰,是最接近随机对照试验的一种研究方法,具有节约时间和成本的特点。其次,MR研究还能避免观察性研究的一些局限性,包括反向因果关系和潜在的混杂因素等。此外,本研究使用了大样本量的纳入研究及与RA显著相关的工具变量,加之多重敏感性分析,从而增强了结果的稳健性。然而研究中也存在一些不足之处,由于研究是基于汇总数据水平的MR分析假设,所以无法研究RA对OP潜在的非线性影响,有必要进一步研究以验证这些关系的性质。除此之外,本研究也未能根据不同的RA和OP类型进行亚组分析,例如原发性OP或继发性OP。最后,考虑到在异质种族人群中等位基因频率可能会导致遗传研究的偏差,所以将MR分析限制在了同质种族人群中,以避免来自更多异质人群的潜在混淆,但对于其他种族人群的解释和应用是有限的。

综上,本研究通过采用孟德尔随机化方法,从遗传学角度探讨了RA与OP的因果关系,为RA与OP之间存在因果关联提供了证据,进而用以指导临床。因此,对于RA患者,可以考虑在疾病早期积极预防和临床干预OP。