骨質疏松癥與潰瘍性結腸炎共病聯系的生物信息學分析

李紹爍 陳柏行 陳浩 尹恒 王建偉*

1.無錫市中醫醫院，南京中醫藥大學無錫附屬醫院，江蘇 無錫 214071

2.比利時魯汶大學骨發育與再生實驗室，比利時，3000

骨質疏松癥(osteoporosis，OP)是一種常見的全身性、代謝性骨骼疾病，其主要特點是骨微結構破壞、骨量下降、骨質脆性增加，導致骨折風險增加[1]。流行病學研究顯示，我國OP患者已逾7 000萬，預計到2050年因骨質疏松癥及骨質疏松性骨折耗費的財政支出將達約250億元[2-3]，骨質疏松癥已成為我國重點攻關研究的老年疾病之一[4]。潰瘍性結腸炎(ulcerative colitis，UC)是一種慢性非特異性炎癥性腸道疾病，是常見的消化系統疑難病，目前UC的發病機制尚不明確，可能與遺傳、免疫、精神心理等因素相關，據推測UC在中國的發病率約為11.6/10萬人[5]。

雖然骨質疏松癥與潰瘍性結腸炎是兩個不同系統的疾病，但人體是一個復雜而統一的整體，臨床證據表明，OP與UC之間存在重要的聯系。骨質疏松癥發病基礎是成骨細胞-破骨細胞調控的骨穩態平衡被打破，成骨細胞調控的骨形成減弱，破骨細胞調控的骨吸收增強。機體炎癥因子異常表達與骨穩態平衡密切相關，研究表明UC患者血清腫瘤壞死因子(TNF)、白介素6(IL-6)、白介素8(IL-8)等炎癥因子水平明顯高于健康人群，而骨密度則明顯低于健康人群[6]，同時，UC病理過程中小腸對維生素D、維生素K、鈣質吸收下降，治療過程中使用糖皮質激素，都可能導致OP的發生或加重[7-8]，而OP發病過程中骨代謝紊亂相關的炎癥因子異常表達、免疫功能失調可能加重UC的病程。

MicroRNA(miRNA)是一類由長度為20～24個核苷酸組成的內源性單鏈非編碼RNA分子，miRNA通過與信使RNA(mRNA)3′非翻譯區(UTR)中的互補序列靶向結合，阻斷蛋白質翻譯，調控轉錄后基因表達，發揮對基因表達的調控作用[9-11]。研究表明，miRNA在調控骨質疏松癥及潰瘍性結腸炎的發病、轉歸、治療中扮演著重要的角色。Li等[12]發現miR-188是調控骨髓間充質干細胞(BMSCs)成骨分化與成脂分化的關鍵位點，可能是老年骨質疏松癥的潛在治療靶點；Suarjana等[13]通過對比絕經后骨質疏松癥婦女與絕經后非骨質疏松癥婦女的血清樣本，發現miR-21在絕經后骨質疏松癥組明顯高表達，同時miR-21與骨保護素(OPG)、核因子κB受體活化因子配體(RANKL)的表達密切相關，認為miR-21激活了破骨細胞的增殖分化活動；Cai等[14]對UC患者進行了miRNA芯片測序，發現UC患者的miR-141可能通過下調趨化因子5(CXCL5)的表達，激活絲氨酸/蘇氨酸激酶(AKT)信號通路，導致UC的發病。隨著現代生物信息技術的高速發展，人類逐漸深入從微觀分子層面認識疾病的發生發展及不同疾病之間的聯系，但目前從“共病”角度探索骨質疏松癥與潰瘍性結腸炎的聯系的研究仍較少。本研究借助現代生物信息學理論方法，通過對比OP與UC的miRNA測序芯片數據，挖掘聯系OP與UC兩個不同疾病的miRNA，嘗試從分子機制層面闡釋OP與UC之間的共病機制，以期為后續研究提供思路方向。

1 資料與方法

1.1 數據獲取

從美國國立生物信息技術中心NCBI的基因表達數據庫GEO(Gene Expression Omnibus，https://www.ncbi.nlm.nih.gov/gds)[15]檢索下載骨質疏松癥患者miRNA數據芯片GSE93883與潰瘍性結腸炎患者miRNA數據芯片GSE63806的原始數據。GSE93883芯片包含12例骨質疏松癥患者與6名健康對照個體的miRNA表達數據，GSE63806芯片包含9例潰瘍性結腸炎患者與12名健康對照個體的miRNA表達數據。

1.2 差異表達miRNA篩選

GEO數據庫在線分析工具GEO2R集成基于R語言的基因表達量數據統計分析功能，可對疾病組與對照組數據進行組間t檢驗統計分析，篩選差異表達基因[16]。本研究使用GEO2R分別對骨質疏松癥與潰瘍性結腸炎miRNA芯片數據進行差異表達分析，設置檢驗統計P<0.05、差異倍數(fold change，FC)的對數絕對值|logFC|>1為篩選條件，分別篩選出OP與UC的差異表達miRNA(differentially expressed miRNAs，DEmiRNAs)，并對DEmiRNAs取交集，獲得聯系OP與UC兩個疾病的關鍵miRNA。

1.3 DEmiRNA靶標基因預測

使用在線miRNA靶標基因預測工具miRWalk2.0(http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2)[17-18]進行OP與UC的共病DEmiRNAs的靶標基因預測。miRWalk集成了包括miRanda[19]、RNA22[20]、miRDB[21]、Targetscan[22]等多個不同的miRNA靶基因預測工具在內，可進行多重數據庫的miRNA共篩選，并找到其中共有的靶標基因，最大限度增加預測置信度。上傳DEmiRNAs至miRWalk2.0數據庫，勾選miRWalk、miRanda、RNA22、miRDB、Targetscan等5個數據庫進行DEmiRNAs的靶標基因預測功能，下載保存預測結果，篩選出5個數據庫共有的靶標基因。

1.4 靶標基因蛋白互作分析

人體生物學功能由復雜的基因表達調控共同完成，基因與基因、非編碼RNA與基因之間存在蛋白產物互作聯系(protein-protein interaction，PPI)，是基因間互相作用的物質基礎。本研究使用在線蛋白互作分析工具STRING V11(https://string-db.org/)[23]進行OP與UC共病miRNA的靶標基因蛋白互作分析，限定物種為“Homo sapiens”(人類)，關聯置信度為<0.40，運行分析，并將結果導入至生物信息學網絡構建分析工具Cytoscape 3.6.1軟件[24-25]進行miRNA靶標基因蛋白互作網絡的構建分析。

1.5 靶標基因生物信息學富集分析

基因本體論(gene ontology，GO)通過對基因產物及功能進行注釋分析了解基因潛在生物學功能[26]，京都基因與基因組百科全書(Kyoto Encyclopedia of Genes and Genomes，KEGG)從分子信號通路水平對基因參與的生物學功能調控作用進行分析[27]，兩者為目前最常見的生物信息學分析方法。本研究使用在線基因生物信息學富集分析數據庫KOBAS3.0(http://kobas.cbi.pku.edu.cn/kobas3)[28]進行靶標基因的GO與KEGG富集分析，設置檢驗統計P<0.05為差異具有統計學意義。

2 結果

2.1 共病miRNAs鑒定

使用GEO2R對骨質疏松癥miRNA芯片GSE93883與潰瘍性結腸炎miRNA芯片GSE63806進行差異表達分析，從GSE93883芯片共篩選獲得424個DEmiRNAs，從GSE63806共篩選獲得11個DEmiRNAs，對兩個芯片的DEmiRNAs取交集，最終獲得1個骨質疏松癥與潰瘍性結腸炎的共病miRNA，為hsa-miR-455-3p。hsa-miR-455-3p在OP與UC芯片的表達情況見表1，取交集情況以韋恩圖(Venn diagram)呈現，見圖1A。

表1 miRNA表達情況Table 1 Expression of miRNA

2.2 miRNA-靶標基因調控網絡構建

研究使用miRWalk2.0集成的多個miRNA靶標基因預測數據庫對hsa-miR-455-3p進行靶標基因預測與收集。從miRDB預測獲得hsa-miR-455-3p靶向于188個基因，miRanda預測獲得hsa-miR-455-3p靶標基因2 057個，Targetscan預測靶標基因3 090個，RNA22預測得靶標基因3 445個，miRWalk預測靶標基因3 142個，將上述5個數據庫獲得的靶基因取共交集，最終獲得hsa-miR-455-3p的靶標基因26個，取共交集情況見圖1B。

圖1 韋恩圖A：OP芯片GSE93883與UC芯片GSE63806取交集；B：miRWalk、miRDB、miRanda、RNA22、Targetscan預測hsa-miR-455-3p靶標基因共交集Fig.1 Venn diagramA: The OP chip gse93883 intersects with UC chip gse63806; B: miRWalk, miRDB, miRanda, RNA22, and Targetscan predict the co-intersection of hsa-mir-455-3p target genes

將獲得的骨質疏松癥與潰瘍性結腸炎共病miRNA(hsa-miR-455-3p)及其預測的26個靶標基因共同導入至Cytoscape軟件，根據其靶向關系構建miRNA-靶標基因調控關系網絡，將miRNA-靶標基因網絡與miRWalk篩選靶標基因結果共同輸出，見圖2。

圖2 靶標基因預測與調控網絡圖Fig.2 Potential targets prediction and miRNA-targets network

2.3 靶標基因蛋白互作網絡

將獲得的hsa-miR-455-3p靶標基因上傳至STRING V11數據庫，運行靶標基因蛋白互作分析，將分析結果下載導入至Cytoscape軟件進行蛋白互作網絡構建，靶標基因蛋白互作網絡見圖3 A，節點(node)分別代表hsa-miR-455-3p的靶標基因，不同節點之間以邊線(edge)相連代表其互作關系。進一步使用Cytoscape軟件的網絡分析工具MCODE對網絡關聯情況進行分析，MCODE共識別出靶標基因蛋白互作網絡中存在3個子網絡，子網絡1(11個節點，55條邊線)、子網絡2(7個節點，18條邊線)、子網絡3(4個節點，5條邊線)，子網絡見圖3B、3C、3D。

圖3 靶標基因互作網絡A：靶標基因互作網絡；B：子網絡1；C：子網絡2；D：子網絡3Fig.3 Protein-protein interaction networkA: Target gene interaction network; B: Sub network 1; C: Sub network 2; D: Sub network 3

2.4 靶標基因GO與KEGG富集分析

將篩選獲得的hsa-miR-455-3p靶標基因整理上傳至KOBAS3.0數據庫，選定物種為“Homo sapiens”(人類)，分析項目為Gene Oncology、KEGG Pathway，運行富集分析并以檢驗統計P<0.05篩選具有統計學差異的結果。結果顯示，26個靶標基因主要顯著富集于8個GO生物學注釋過程，18個KEGG信號通路，將GO與KEGG富集分析結果整理匯總，見表2、3及圖4。

表2 GO富集分析結果Table 2 Results of GO enrichment

表3 KEGG信號通路富集分析結果Table 3 Results of KEGG pathways enrichment

圖4 GO與KEGG富集分析結果A：GO富集分析結果，顏色越深代表富集程度越高，氣泡大小代表富集基因數目；B：KEGG信號通路富集分析結果，顏色越深代表富集程度越高，氣泡大小代表富集基因數目Fig.4 Results of GO enrichment and KEGG pathways enrichmentA: The deeper the color, the higher the enrichment degree, the bubble size, and the number of enrichment genes; B: The results of KEGG signal pathway enrichment analysis showed that the darker the color, the higher the enrichment degree, the bubble size, and the number of enriched genes

3 討論

探索骨質疏松癥與潰瘍性結腸炎發生、發展背后的分子機制將有助于疾病后續的診斷、治療及預后。基因芯片等高通量測序技術的飛速發展，為我們從微觀分子水平認知疾病提供了機會，更為尋找疾病新的治療靶點提供了客觀基礎。研究表明，miRNA參與調控約占人體60%的基因組轉錄翻譯，廣泛參與目前認知的細胞生物學過程的調節，其表達水平與許多疾病的發生發展密切相關[29]。本研究立足于對骨質疏松癥與潰瘍性結腸炎臨床證據的聯系，運用現代生物信息學技術方法，對GEO數據庫中收錄的OP和UC的miRNA測序芯片數據進行挖掘分析，獲得OP與UC的共病聯系DEmiRNA，并借助在線預測工具獲得miRNA調控的靶標基因，對靶標基因進行下游相關的生物學過程、信號通路、基因蛋白互作等分析，探索OP與UC之間潛在的分子機制聯系。

借助對OP基因芯片GSE93883與UC基因芯片GSE63806進行差異分析及交集映射，獲得了聯系兩疾病的共病DEmiRNA——hsa-miR-455-3p，目前探索has-miR-455-3p在人體生命活動中所起的功能作用展開的研究已有不少。研究發現hsa-miR-455-3p在骨代謝活動、炎癥反應、細胞周期調控中扮演重要角色，Zhang等[30]使用枸櫞酸鐵銨(FAC)干預成骨細胞造成細胞內鐵過載，建立骨質疏松癥模型，發現hsa-miR-455-3p可以通過下調HDAC2蛋白表達，調控Nrf2/ARE信號通路，抑制氧化應激反應，促進成骨細胞增殖；Zeng等[31]發現hsa-miR-455-3p可抑制TGF-β/Smad信號通路關鍵基因ZEB1，調控TGF-β/Smad信號通路，維持細胞分化活動；Hu等[32]發現hsa-miR-455-3p可靶向于TGF-β/Smad信號通路的PAK2基因，抑制軟骨退變；Zhu等[33]發現hsa-miR-455-3p調控Wnt/β-Catenin信號通路抑制高糖誘導的系膜細胞炎癥反應；Fiorillo等[34]通過體內及體外實驗驗證，證明hsa-miR-455-3p參與調控NF-κB信號通路，發揮類似于糖皮質激素的抑制慢性炎癥的作用。

miRNA廣泛參與人類、動植物基因的轉錄后調控，以此參與細胞增殖、分化、凋亡等復雜生命活動進程，在這個過程中，miRNA與其靶標基因的結合是發揮其調控作用的物質基礎。本研究中通過使用miRWalk、miRDB、RNA22、miRanda、Targetscan等5個miRNA靶標基因預測數據庫，獲取了hsa-miR-455-3p的26個靶標基因，并構建了miRNA-靶標基因調控網絡，進一步分析發現靶標基因間存在復雜網絡調控的蛋白互作關系。hsa-miR-455-3p的多個靶標基因目前已被證實通過多種信號通路活動參與OP與UC的病理活動。靶標基因p21活化蛋白激酶2(PAK2)被認為是咖啡因代謝導致的骨丟失中的重要靶點，Lu等[35]通過建立咖啡因誘導成骨細胞凋亡的實驗，觀察發現PAk2是其中的中介位點，咖啡因代謝激活PAK2促進成骨細胞凋亡；TGF-β/Smad信號通路是經典的骨代謝調控通路，靶標基因Smad2在轉導細胞膜表面受體的TGF-β信號進入胞核中發揮重要作用，Qiu等[36]觀察發現miRNA-214可調控TGF-β/Smad2促進骨髓間充質干細胞成脂分化，抑制成骨分化活動；靶標基因XDH編碼黃嘌呤脫氫酶，黃嘌呤脫氫酶是參與氧化代謝的重要羥化酶成員，一項關于炎癥性腸病的基因篩查研究發現，XDH基因多態性變化影響了UC患者對免疫抑制治療的敏感性[37]。

同時，GO、KEGG富集分析結果顯示，hsa-miR-455-3p的靶標基因主要富集于細胞代謝活動、細胞周期、咖啡因代謝等生物學過程及信號通路。細胞進程(cellular process)、細胞成分組成(cellular component organization)主要指微觀細胞層面的結構特點，含堿基化合物代謝過程正向調控(positive regulation of nucleobase-containing compound metabolic process)、離子跨膜運輸(anion transmembrane transport)等主要指細胞的功能活動，均為細胞基礎生命活動形式，hsa-miR-455-3p可能通過調控靶標基因在上述細胞基礎活動中的表達，影響疾病的病理機制。叉行頭轉錄因子FoxO信號通路是控制細胞凋亡的重要通路，其上游受到PI3K/AKT信號通路、NF-κB信號通路等多重調控，PI3K/AKT信號通路的激活可上調成骨細胞ALP、BMP2等表達，促進成骨細胞增殖分化，同時還可抑制局部炎癥活動[38-39]；嘌呤代謝、咖啡因代謝等活動與骨代謝及骨細胞增殖分化密切相關，嘌呤、咖啡因等物質也與氧化應激反應、炎癥因子聚集反應直接相關，高尿酸、高咖啡因攝入被認為是骨質疏松癥的危險因素[40]。

本研究中篩選獲得的OP與UC共病miRNA(hsa-miR-455-3p)在骨質疏松癥患者中表達下調，但在潰瘍性結腸炎患者中表達上調，呈現相反的表達趨勢，hsa-miR-455-3p的靶標基因在調控OP與UC相關的骨代謝、炎癥反應、細胞周期等過程中亦表現出或相同或不同的趨勢，提示OP與UC之間存在深層次復雜聯系，并非簡單的正相關或負相關調控關系。OP與UC均為病理機制尚不完全明確的復雜疾病，與代謝活動、炎癥反應、免疫調控等多重因素密切相關，本研究以miRNA為聯系OP與UC的切入點，通過挖掘OP與UC患者的miRNA基因芯片，篩選出對OP與UC存在共調控作用的差異miRNA，發現該miRNA調控的靶標基因主要通過參與成骨細胞增殖分化活動、細胞轉運活動、嘌呤及咖啡因代謝等環節影響OP與UC的疾病發生發展進程，認為hsa-miR-455-3p可能是聯系OP與UC的關鍵點，是治療OP與UC的潛在治療靶點。另一方面，本研究僅從個別OP與UC的基因芯片數據中篩選獲得聯系兩疾病的miRNA，缺少更大數據樣本量的研究支持。目前關于miRNA的研究仍處于起步階段，存在許多未知的miRNA及其未知的功能等待人類挖掘發現，希望今后能繼續深入探索miRNA在不同疾病間的聯系，幫助從整體對不同疾病病理機制深入理解，服務于臨床。