基于孟德爾隨機化研究類風濕性關節炎與骨質疏松癥的因果關系

朱凱 劉宛欣 王秋根

上海中醫藥大學附屬岳陽中西醫結合醫院,上海 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。