基于生物信息學(xué)探討骨質(zhì)疏松癥與膝骨關(guān)節(jié)炎的關(guān)系

陳桐瑩 高豐禾 汪悅東 林燕平 黃佳純 黃宏興 萬雷* 王世浩 趙宇 劉遠(yuǎn)祥 袁嘉堯

1.廣州中醫(yī)藥大學(xué)第三臨床醫(yī)學(xué)院，廣東 廣州 510405 2.廣州中醫(yī)藥大學(xué)第三附屬醫(yī)院，廣東 廣州 510375

骨質(zhì)疏松癥(osteoporosis，OP)是常見的骨骼疾病，主要特征是骨丟失增加和骨微結(jié)構(gòu)惡化，導(dǎo)致骨強度降低和骨折風(fēng)險增加。早期流行病學(xué)調(diào)查顯示我國OP患者近7 000萬，骨量減少者已超過2億人[1]。2018年的OP流行病學(xué)調(diào)查顯示，到2050年，OP性骨折數(shù)量將增加到599萬例，耗費254.3億元[2]。膝骨關(guān)節(jié)炎(knee osteoarthritis，KOA)是一種緩慢進行性的變性關(guān)節(jié)疾病，主要特征是關(guān)節(jié)疼痛和關(guān)節(jié)軟骨的進行性退行性變。據(jù)統(tǒng)計，中國大陸地區(qū)患KOA人數(shù)已超過1億[3]。OP和KOA均是影響人們生活質(zhì)量的慢性骨骼肌肉疾病，會為個人和國家?guī)沓林氐呢?fù)擔(dān)。

OP和KOA本是兩種相互獨立的疾病，但臨床表明OP和KOA可在同一個體共存并互相影響，加劇兩病的進展[4]。關(guān)于兩者的關(guān)系，基礎(chǔ)與臨床研究數(shù)量逐漸增加。目前研究主要集中在軟骨下骨特性、遺傳、BMI與BMD之間相互作用方面，其中一些因素對兩病影響是相同的(如年齡、性別、遺傳學(xué)和炎癥)，而其他的則是相反(BMI、BMD)[5]，但OP和KOA的關(guān)系仍然不明確，仍需要更深入地剖析。目前發(fā)現(xiàn)OP和KOA間基因易感性有一定的聯(lián)系，如α-2-HS糖蛋白、M型肌酸激酶在OP與KOA中均存在密切聯(lián)系[6]。這提示可以通過兩病中高敏感性與高特異性的交集基因探索OP和KOA的關(guān)系。

MicroRNA(miRNA)是短內(nèi)源性非編碼RNA分子，長度為19～25個核苷酸，可通過直接靶向基因的3'-UTR來調(diào)控轉(zhuǎn)錄后基因的表達(dá)，對基因表達(dá)產(chǎn)生負(fù)調(diào)控[7]。研究表明miRNA在OP和KOA的發(fā)生發(fā)展中起著重要的作用。如miR-141抑制Calcr和EphA2的表達(dá)進而減輕恒河猴的骨吸收減輕OP癥狀[8]，miR-103通過靶向Sox6促進OA的發(fā)展[9]，miR-29靶向抑制VEGF減輕OA癥狀等[10]。隨著生物信息的快速發(fā)展，諸多疾病和基因的數(shù)據(jù)庫逐漸完善，目前探索OP和KOA關(guān)系的研究很少深入到非編碼RNA方面，且現(xiàn)有關(guān)于OP和KOA miRNA生物信息學(xué)研究和對miRNA靶向核心基因報道更為稀缺。

本研究通過對基因芯片GSE93883、GSE105027進行生物信息學(xué)分析篩選核心交集差異miRNA，并預(yù)測與這些核心miRNA的靶向基因，為闡明OP和KOA的關(guān)系以及后續(xù)更深入的研究提供依據(jù)，并為OP和KOA治療藥物的研發(fā)提供較為可靠的路徑和作用靶點。

1 材料與方法

1.1 數(shù)據(jù)收集

從GEO(Gene Expression Omnibus database)公共數(shù)據(jù)庫(http://www.ncbi.nlm.nih.gov/geo)下載骨質(zhì)疏松癥血清基因表達(dá)芯片GSE93883，包含12位骨質(zhì)疏松患者和6位健康個體，其所處平臺為GPL21575 Agilent-070156 Human miRNA V21.0 Microarray 046064；膝骨關(guān)節(jié)炎血清基因表達(dá)芯片為GSE105027，共含12位原發(fā)性KOA患者和12位健康個體，其所處平臺為GPL21575 Agilent-070156 Human miRNA V21.0 Microarray 046064。

1.2 差異miRNA的數(shù)據(jù)分析

使用GEO2R在線分析軟件進行差異miRNA (differentially expressed miRNAs，DEmiRNAs)的篩選[11]，GEO2R是GEO數(shù)據(jù)庫自帶的在線分析工具，通過R語言分析數(shù)據(jù)。經(jīng)原始數(shù)據(jù)進行成組t檢驗統(tǒng)計學(xué)分析，以adj.P<0.05和|logFC|≥1作為篩選條件。使用Venny 2.1篩選出OP和KOA交集DEmiRNAs。

1.3 靶基因預(yù)測

使用TargetScan(http://www.targetscan.org/vert_71/)和miRDB(http://mirdb.org)數(shù)據(jù)庫分別預(yù)測了DEmiRNA的靶基因，使用Venny生成交集基因。

1.4 基因富集分析

GO(Gene Ontology)是基因分析的常用方法，可以對基因進行注釋和分類。GO分析由分子功能(molecular function)、生物過程(biological process)和細(xì)胞組成(cellular component)3個部分組成。KEGG(kyoto encyclopedia of genes and genomes)通路也是常見的基因功能富集分析方法。本研究應(yīng)用DAVID數(shù)據(jù)庫(Database for Annotation，Visualization and Integrated Discovery)進行在線分析提供所需的 GO和KEGG生物功能富集數(shù)據(jù)[12]。本研究使用的DAVID數(shù)據(jù)庫版本為6.8，地址為https://david.ncifcrf.gov/。GO各項、KEGG各項均以P<0.05為篩選條件。

1.5 蛋白互作網(wǎng)絡(luò)與模型分析

Cytoscape(https://cytoscape.org/)主要用于整合模塊化網(wǎng)絡(luò)和生物科學(xué)聯(lián)系網(wǎng)絡(luò)圖的繪制。本研究將篩選的DEmiRNAs與篩選出的靶向基因植入cytoscape(版本3.7.1)，以構(gòu)建miRNA-mRNA網(wǎng)絡(luò)。STRING(search tool for the retrieval of interacting genes)(https://string-db.org/)是用于評估和制作蛋白互作網(wǎng)絡(luò)的在線分析工具[13]。本研究將篩選到的靶向基因?qū)隨TRING數(shù)據(jù)庫探索基因間的潛在聯(lián)系。此后，把STRING的計算結(jié)果導(dǎo)入Cytoscape進行MCODE分析以挖掘PPI中連接最為緊密的集簇。使用cytohubba插件進行計算，最終通過DEGREE法篩選出top10的基因。

2 結(jié)果

2.1 DEmiRNA的鑒定

本研究選用基因表達(dá)芯片GSE93883、GPL21575。經(jīng)GEO2R初步分析獲得54 107個DEmiRNAs，隨后由經(jīng)P<0.05，|logFC|≥1條件篩選并取交集，最終獲4個DEmiRNAs：hsa-miR-3202、hsa-miR-4687-3p、hsa-miR-4508、hsa-miR-320b。其中OP組均上調(diào)，而OA組均下調(diào)。見表1。

表1 交集DEmiRNAsTable 1 Differential expression of miRNAs

2.2 miRNA-基因調(diào)控網(wǎng)絡(luò)預(yù)測

本研究運用TargetScan和miRDB數(shù)據(jù)庫預(yù)測上述篩選的4個DEmiRNAs靶向基因，結(jié)果顯示共有12 141個靶向基因，Venny分別篩選靶向基因在兩個數(shù)據(jù)庫的交集，結(jié)果發(fā)現(xiàn)hsa-miR-3202靶向248個基因，hsa-miR-4687-3p靶向79個基因，hsa-miR-4508靶向5個基因，hsa-miR-320b靶向25個基因。在cytoscape構(gòu)建miRNA-基因調(diào)控網(wǎng)絡(luò)圖，圖形越大，顏色越深，Degree得分越高，見圖1。

圖1 交集miRNA-mRNA調(diào)控網(wǎng)絡(luò)Fig.1 miRNA-mRNA regulatory network

2.3 GO功能與KEGG信號通路富集分析

把篩選的基因?qū)隓AVID數(shù)據(jù)庫，僅選用“Homo sapiens”，分別分析獲得疾病基因和蛋白信息所富集的通路，并將基因信息可視化。在GO分析上，本研究以P<0.05為篩選條件。結(jié)果發(fā)現(xiàn)在生物學(xué)過程中，靶基因主要聚集在轉(zhuǎn)錄、細(xì)胞內(nèi)蛋白質(zhì)運輸、RNA聚合酶正負(fù)調(diào)控上(表2)；在細(xì)胞定位中，基因主要富集在細(xì)胞質(zhì)和核質(zhì)中(表3)。從分子功能上看，基因主要集中在蛋白質(zhì)結(jié)合和核酸結(jié)合上(表4)。KEGG分析中本研究以P<0.05為篩選條件，結(jié)果顯示基因主要集中在癌癥通路、Rap1信號通路、鈣信號通路、Pi3k-Akt、Wnt通路、cAMP通路等(表5)。

表2 差異靶基因生物學(xué)過程分析Table 2 Biological process analysis of the target genes

表3 差異靶基因細(xì)胞組件分析Table 3 Cellular component analysis of the target genes

2.4 核心基因的鑒定與篩選的集簇模塊

通過 STRING的PPI構(gòu)建(圖2)，經(jīng)由 Cytoscape對網(wǎng)絡(luò)的計算工具得出所有靶向基因的連接度(degree)。本研究按degree從高至低進行排序，以排名前10位的高degree 基因定為核心基因，它們分別是VEGFA(degree=30)、ESR1(degree=28)、CCND1(degree=27)、PRKACA(degree=23)、GSK3β(degree=21)、H2AFX(degree=17)、SHH(degree=17)、EGR1(degree=16)、DNMT3A(degree=15)、DNMT3B(degree=14)(圖3，圖中紅色代表cytohubba插件得分最高，黃色次之)。此外，MCODE分析發(fā)現(xiàn)11組集簇，共包含61個節(jié)點(node)與127條連線(edge)。本研究以score為依據(jù)展示前 3組集簇并分析各個集簇所富集的通路(圖4)。

表4 差異靶基因分子功能分析Table 4 Functional analysis of the target genes

表5 靶向基因的KEGG富集分析Table 5 KEGG pathway analysis of the target genes

圖2 靶基因的PPI網(wǎng)絡(luò)Fig.2 The PPI network of the target genes

圖3 前10hub基因Fig.3 Top 10 hub genes

3 討論

圖4 蛋白質(zhì)互作網(wǎng)絡(luò)的前3個集簇模塊Fig.4 Top three modules from the PPI network

MiRNA作為一類基因表達(dá)調(diào)節(jié)劑，在正常的骨組織生理及病理學(xué)中都非常重要。研究表明miR-138、miR-338-3p和miR-188均可以抑制OP的發(fā)生[14]。miR-125b通過靶向BMPR2調(diào)節(jié)軟骨細(xì)胞的增殖分化[15]。本研究通過對GEO數(shù)據(jù)庫中的基因芯片GSE93883和GSE105027進行生物信息學(xué)分析獲得DEmiRNAS，預(yù)測靶基因并分析了有關(guān)這些基因富集的生物過程、細(xì)胞定位、分子功能和信號通路。其中hsa-miR-3202靶向的基因為VEGFA、CCND1、GSK3β、SHH、DNMT3A、DNMT3B。hsa-miR-4687-3p則靶向ESR1、EGR1，以上均為篩選得到的前十基因，這說明hsa-miR-3202和hsa-miR-4687-3p與OP和KOA的發(fā)病關(guān)系較緊密，可以成為后續(xù)研究的對象，為闡明OP和KOA的關(guān)系以及后續(xù)更深入的研究提供依據(jù)。

MiRNA是一個龐大復(fù)雜的網(wǎng)絡(luò)調(diào)控體系，對miRNA相關(guān)通路的研究有助于了解其調(diào)控機理的重要方式。GO與 KEGG分析有助于更深入地認(rèn)識并篩選出DEmiRNA靶向基因的功能和作用。本研究由KEGG篩選出的富集通路以癌癥通路為首，其次為PI3K-Akt、鈣信號通路、Rap1信號通路、Wnt通路、cAMP通路等。其中PI3K-Akt通路與OP及KOA有著密切的聯(lián)系，PI3K/AKT信號的激活轉(zhuǎn)導(dǎo)上調(diào)成骨細(xì)胞分化標(biāo)記基因的表達(dá)，包括ALP和BMP-2，從而促進成骨細(xì)胞的增殖和分化[16]。研究發(fā)現(xiàn)PI3K/AKT/mTOR途徑在延遲炎癥中起關(guān)鍵作用，并且還調(diào)節(jié)軟骨細(xì)胞凋亡、代謝、基因表達(dá)等[17]，阻斷PI3K/AKT信號傳導(dǎo)可抑制軟骨下骨的骨硬化并減輕創(chuàng)傷后OA的發(fā)生[18]。Rap1是一種小的GTP酶，通過調(diào)節(jié)細(xì)胞間基質(zhì)粘附和肌動蛋白重排，動態(tài)控制細(xì)胞擴散和內(nèi)皮屏障功能的協(xié)調(diào)[19]。Rap1/MAPK信號通路可促進成骨細(xì)胞增殖分化并抑制凋亡[20]。而在OA此通路研究較少，可成為探索兩病關(guān)系的切入點。Wnt信號被認(rèn)為是骨骼生物學(xué)中最重要的途徑，在骨形成中起著關(guān)鍵作用，而Wnt信號轉(zhuǎn)導(dǎo)的上調(diào)卻可能是KOA發(fā)展和進程的一個致病因素[21]。兩者之間的矛盾關(guān)系也是探索兩病關(guān)系的切入點之一。

MiRNA功能研究的重點在于對其靶基因調(diào)控的研究。本研究最終篩選出10個miRNA的靶基因：VEGFA、ESR1、CCND1、PRKACA、GSK3Β、H2AFX、SHH、EGR1、DNMT3A、DNMT3B。其中VEGFA即血管內(nèi)皮生長因子，是血管生成和成骨的有效耦合所必需的。VEGF可以部分通過自分泌和內(nèi)分泌機制以及旁分泌信號調(diào)節(jié)成骨細(xì)胞、破骨細(xì)胞的功能，起到促進骨形成減少骨吸收的作用[22]。而OA患者中VEGF表達(dá)水平上升，OA隨即加重[23]。ESR1是雌激素受體，GWAS揭示了ESR1是與骨密度(BMD)相關(guān)的新型遺傳基因座，其能影響骨質(zhì)疏松的發(fā)展[24]。由于雌激素可能同時具有抗炎和促炎特性，所以ESR1是否會加劇骨關(guān)節(jié)炎的發(fā)生仍有爭議[25]。CCND1是細(xì)胞周期蛋白D1，主要功能是調(diào)節(jié)細(xì)胞周期從DNA的預(yù)合成階段到合成階段。CCND1的異常表達(dá)極大地影響細(xì)胞周期進程。研究發(fā)現(xiàn)綠原酸通過PI3K/Akt/CCND1途徑預(yù)防去卵巢大鼠的骨質(zhì)疏松癥[26]。CCND1基因沉默則會抑制軟骨細(xì)胞增殖[27]。GSK3β是Wnt通路的關(guān)鍵蛋白，GSK3β降解，則β-catenin磷酸化降低繼而促進成骨作用[28]。而GSK3β活性在OA中卻具有軟骨保護作用[29]。SHH是Hedgehog的同源基因，研究表明SHH能促成骨分化[30]，相反地SHH加劇OA中的軟骨降解[31]。EGR1是表皮生長因子轉(zhuǎn)錄因子之一，可促進成骨細(xì)胞增殖，減少細(xì)胞凋亡[32]。而研究發(fā)現(xiàn)EGR1降低可能提示OA處于早期狀態(tài)[33]。DNMT是DNA甲基轉(zhuǎn)移酶的一種，可引起DNA甲基化進而抑制靶基因，有研究在缺氧和氧化應(yīng)激誘導(dǎo)的染色質(zhì)修飾酶表達(dá)變化的體外觀察中發(fā)現(xiàn)，DNMT3A的表達(dá)在OP和OA中均顯著下調(diào)[34]。另有研究發(fā)現(xiàn)DNMT3B是軟骨細(xì)胞成熟和軟骨內(nèi)骨形成所必需的[35]，而在OP中DNMT3B卻起著抑制成骨細(xì)胞增殖的作用[36]。

本研究發(fā)現(xiàn)OP和KOA交集的miRNA趨勢相反，且預(yù)測靶向基因中大部分對兩病的作用機制也是截然相反，這說明可能在血液檢測上，骨質(zhì)疏松癥與膝骨關(guān)節(jié)炎存在負(fù)相關(guān)的關(guān)系，但由于本研究選取的數(shù)據(jù)有一定的局限，且骨質(zhì)疏松癥和膝骨關(guān)節(jié)炎均屬復(fù)雜的多因素的代謝性骨骼肌肉疾病，例如Gun認(rèn)為KOA患者較高的BMD并不能證明KOA與OP呈負(fù)相關(guān)關(guān)系，因為骨關(guān)節(jié)炎引起的姿勢不穩(wěn)和肌肉衰弱可以抵消骨關(guān)節(jié)炎患者骨中小梁分離減少和骨小梁數(shù)目增加的結(jié)構(gòu)優(yōu)勢，在此狀態(tài)下，KOA很可能加劇OP的發(fā)生[37]。由此可見，兩病的關(guān)系要從多方面探索并綜合評價。同時，本研究結(jié)果表明OP和KOA交集miRNA可能通過調(diào)控PI3K-Akt、鈣信號通路、Rap1信號通路、Wnt通路、cAMP通路對兩病起作用，這為后續(xù)探索OP與KOA的關(guān)系提供了更多的思路。本研究通過對差異 miRNA 及其調(diào)控的靶基因深入研究，將有助于探索OP和KOA的關(guān)系及為治療提供新的思路與方法，并為OP和KOA治療藥物的研發(fā)提供較為可靠的路徑和作用靶點。