柴毅 譚峰 樊巧玲

南京中醫(yī)藥大學(xué)，江蘇 南京 210046

骨質(zhì)疏松癥(osteoporosis，OP)是最常見(jiàn)的骨骼疾病，以骨量低，骨組織結(jié)構(gòu)破壞，最終導(dǎo)致骨脆性增加，骨強(qiáng)度下降及骨折風(fēng)險(xiǎn)增加，易發(fā)生骨折為特征的全身性骨病[1]。在美國(guó)，每年約有150萬(wàn)人次發(fā)生骨折，絕大多數(shù)都發(fā)生于絕經(jīng)后婦女[2]。絕經(jīng)后骨質(zhì)疏松癥(postmenopausal osteoporosis，PMOP)是與年齡相關(guān)的衰老性疾病，多發(fā)生于絕經(jīng)2年以上，70歲以下的婦女。PMOP的病變是一個(gè)隱性的過(guò)程，全世界約有50%的絕經(jīng)后婦女受到影響，被認(rèn)為是老年人發(fā)病率最高的疾病之一[3]。PMOP已成為全球經(jīng)濟(jì)的負(fù)擔(dān)之一，積極開(kāi)展對(duì)PMOP的預(yù)防和治療是公共衛(wèi)生的重要任務(wù)。

OP的診斷基于全面的病史記錄、體格檢查、骨密度測(cè)定、影像學(xué)檢查和必要的生化測(cè)定，OP的診斷主要基于雙能X線吸收檢測(cè)法(dual-energy X-ray absorptiometry，DXA)骨密度測(cè)量結(jié)果與脆性骨折[4]。然而DXA檢測(cè)也存在一些問(wèn)題，例如DXA檢測(cè)成本較高，不同設(shè)備的DXA檢測(cè)結(jié)果存在差異[5-6]。挖掘與OP或PMOP病理相關(guān)的核心基因，特別是高敏感性與高特異性的基因是預(yù)防和治療該疾病性價(jià)比較高的途徑之一。早期的研究發(fā)現(xiàn)了一些與OP有關(guān)的基因。例如，細(xì)胞周期蛋白E1(Cyclin E1，CCNE1)是細(xì)胞周期的調(diào)控因子。CCNE1參與了骨代謝過(guò)程，CCNE1在PMOP B細(xì)胞中的表達(dá)呈下調(diào)趨勢(shì)[7]。又如細(xì)絲蛋白α(filamin A alpha，F(xiàn)LNA)是參與破骨細(xì)胞生成過(guò)程的關(guān)鍵因子，高表達(dá)的FLNA可促進(jìn)破骨細(xì)胞生成[8]。微小核糖核酸(micro-ribonucleic acid，miRNA)是一類(lèi)非編碼蛋白的小RNA，miRNA通過(guò)抑制特定靶點(diǎn)的mRNA調(diào)控基因表達(dá)[9]。一些miRNA亦可作為OP或PMOP的敏感標(biāo)志物或治療藥物的靶點(diǎn)，因而針對(duì)miRNA在OP或PMOP中的機(jī)制挖掘是很有必要的[10-11]。目前，治療OP或PMOP的相關(guān)機(jī)制的探索很少深入到非編碼RNA的作用層面。此外，利用生物信息學(xué)研究方法在OP或PMOP方面有關(guān)核心基因挖掘的研究數(shù)量相對(duì)較少，且現(xiàn)有關(guān)于鑒定OP或PMOP核心基因的生物信息學(xué)研究對(duì)與核心基因相關(guān)互作的miRNA預(yù)測(cè)的報(bào)道更為稀缺。

本研究通過(guò)對(duì)基因芯片GSE57273進(jìn)行生物信息學(xué)分析篩選核心基因，并預(yù)測(cè)與這些核心基因相互作用的miRNA，為PMOP建立新的科學(xué)假說(shuō)以及后續(xù)更深入的研究提供依據(jù)，并為PMOP診斷以及治療藥物的研發(fā)提供較為可靠的路徑和作用靶點(diǎn)。

1 材料與方法

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

從GEO公共數(shù)據(jù)庫(kù)(https://www.ncbi.nlm.nih.gov/geo/)下載基因表達(dá)芯片GSE57273。該芯片包含3組PMOP藥物干預(yù)前后樣本。其所處的平臺(tái)(Platforms)為GPL4133Agilent-014850 Whole Human Genome Microarray 4x44 K G4112F(Feature Number version)。另下載Series Matrix File以便后續(xù)使用。

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

使用GEO2R(https://www.ncbi.nlm.nih.gov/geo/geo2r/)和Morpheus(https://software.broadinstitute.org/morpheus/)在線分析軟件進(jìn)行差異基因(differentially expressed genes，DEGs)的甄別和篩選。GEO2R是GEO數(shù)據(jù)庫(kù)自帶的公共在線分析工具，它可將GEO數(shù)據(jù)進(jìn)行復(fù)雜的R語(yǔ)言分析，從而呈現(xiàn)出每個(gè)基因的計(jì)算結(jié)果。經(jīng)原始數(shù)據(jù)進(jìn)行成組t檢驗(yàn)統(tǒng)計(jì)學(xué)分析，以adj.P<0.01和|logFC|≥3作為DEGs的篩選條件[12]。

1.3 基因功能富集分析

GO(gene ontology)是常用的分析方法，其主要功能是注釋基因或其產(chǎn)物并識(shí)別高通量基因組或轉(zhuǎn)錄組數(shù)據(jù)的特征生物學(xué)特性。GO按照生物途徑(biological process，BP)、分子功能(molecular function，MF)、細(xì)胞定位(cellular component，CC)對(duì)基因進(jìn)行注釋和分類(lèi)。此外，KEGG(Kyoto Encyclopedia of Genes and Genomes)數(shù)據(jù)庫(kù)可供查詢通路信息和信號(hào)通路檢索等。KEGG通路分析是另一種常用的基因功能富集分析方法。本研究應(yīng)用DAVID(Database for Annotation, Visualization and Integrated Discovery)進(jìn)行在線分析提供所需的GO和KEGG生物功能富集數(shù)據(jù)。本研究使用的DAVID數(shù)據(jù)庫(kù)版本為6.8，地址為https://david.ncifcrf.gov/，由美國(guó)國(guó)立變態(tài)反應(yīng)與傳染病研究所提供研究服務(wù)。使用Fisher Exact或EASE Score統(tǒng)計(jì)方法，GO各項(xiàng)以P<0.05且FDR<0.05為篩選條件，KEGG各項(xiàng)以P<0.05為篩選條件[13]。

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

STRING(Search Tool for the Retrieval of Interacting Genes)是一款可以用來(lái)呈現(xiàn)與評(píng)估蛋白互作(protein-protein interaction，PPI)的在線分析工具。STRING中的所有數(shù)據(jù)和下載文件都可以在“Creative Commons BY 4.0”許可下免費(fèi)獲取。本研究將所篩選的所有DEGs植入STRING(版本10.5，https://string-db.org/)分析工具試探它們之間潛在的聯(lián)系。置信度(confidence score)≥0.4，互作最大值(maximum number of interactors)=0設(shè)為篩選條件[14]。此后，把STRING的計(jì)算結(jié)果導(dǎo)入Cytoscape(版本3.6.0)進(jìn)行MCODE(Molecular Complex Detection)分析以挖掘PPI中連接最為緊密的集簇。本研究使用的MCODE版本為1.5.1，設(shè)置參數(shù)為degree=2，node score=0.2，k-core=2，max. depth=100[15]。

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

CyTargetLinker可以擴(kuò)展生物調(diào)控互作網(wǎng)絡(luò)(regulatory interaction networks，RegINs)，由荷蘭系統(tǒng)生物學(xué)聯(lián)合會(huì)提供支持。它涵蓋了miRNA—靶點(diǎn)、轉(zhuǎn)錄因子—靶點(diǎn)和藥物—靶點(diǎn)之間的互作關(guān)系。本研究下載了人類(lèi)物種基因數(shù)據(jù)集(https://projects.bigcat.unimaas.nl/cytargetlinker/regins/)。選用該數(shù)據(jù)集中基于實(shí)驗(yàn)驗(yàn)證的miRTarBase 4.4數(shù)據(jù)庫(kù)(含20 942個(gè)RegINs)，基于預(yù)測(cè)功能的TargetScan 6.2數(shù)據(jù)庫(kù)(含511 040個(gè)RegINs)和MicroCosm 5數(shù)據(jù)庫(kù)(含541 039個(gè)RegINs)預(yù)測(cè)核心基因與miRNA之間的調(diào)控關(guān)系。

2 結(jié)果

2.1 DEGs的鑒定

本研究選用基因表達(dá)芯片GSE57273，經(jīng)GEO2R初步分析共獲得32 996個(gè)DEGs，隨后由Morpheus分析并經(jīng)條件篩選，最終獲得841個(gè)DEGs，其中包含826個(gè)下調(diào)基因和15個(gè)上調(diào)基因。

2.2 GO功能與KEGG信號(hào)通路富集分析

根據(jù)GO的分析結(jié)果，本研究以P<0.05，F(xiàn)DR<0.05為篩選條件，并按照計(jì)數(shù)值從大到小排列，在BP、CC與MF類(lèi)別中各選取前3項(xiàng)列為表1?？梢钥闯?，在生物學(xué)過(guò)程中，這些DEGs主要參與了基因表達(dá)，細(xì)胞大分子生物合成和RNA代謝過(guò)程；在細(xì)胞定位中，這些DEGs富集于核漿、細(xì)胞質(zhì)基質(zhì)以及粘附連接；從分子功能上看，這些DEGs具有使有機(jī)環(huán)狀化合物結(jié)合、雜環(huán)化合物結(jié)合和核酸結(jié)合的作用。根據(jù)KEGG分析結(jié)果，本研究以P<0.05為篩選條件，按照計(jì)數(shù)值從大到小將DEGs富集的信號(hào)通路列為表2。結(jié)果顯示，這些DEGs主要富集于癌癥信號(hào)通路、病毒致癌通路、粘附斑、rap1信號(hào)通路和內(nèi)質(zhì)網(wǎng)蛋白加工通路。

表1 與PMOP相關(guān)的DEGs的GO富集分析

表2 與PMOP相關(guān)的DEGs的KEGG富集分析

圖2 蛋白質(zhì)互作網(wǎng)絡(luò)的前3個(gè)集簇模塊

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

通過(guò)STRING的PPI構(gòu)建，經(jīng)由Cytoscape對(duì)網(wǎng)絡(luò)的計(jì)算工具得出所有DEGs的連接度(degree)(圖1)。degree值表示網(wǎng)絡(luò)中某一基因與周?chē)虻年P(guān)系數(shù)量，因此degree越大代表與它相互作用關(guān)系的基因數(shù)量就越多。本研究按degree從高至低進(jìn)行排序，以排名前10位的高degree基因定為核心基因，它們分別是HSP90AA1(degree=75)、EP300(degree=55)、SMARCA2(degree=44)、RANBP2(degree=41)、ASH1L(degree=36)、EIF4E(degree=35)、PTEN(degree=31)、CNOT6L(degree=30)、RPL7(degree=29)、KRAS(degree=29)。此外，MCODE分析發(fā)現(xiàn)17組集簇，共包含523個(gè)節(jié)點(diǎn)(node)與2 026條連線(edge)。本研究以score為依據(jù)展示前3組集簇并分析各個(gè)集簇所富集的通路(圖2)。其中集簇模塊A富集核糖體(hsa03010：Ribosome)與mRNA監(jiān)視通路(hsa03015：mRNA surveillance pathway)(P<0.05)，集簇模塊B富集泛素介導(dǎo)的蛋白質(zhì)水解4通路(hsa04120：Ubiquitin mediated proteolysis 4)(P<0.05)，而集簇模塊C未鑒定出具有統(tǒng)計(jì)學(xué)意義的通路(P>0.05)(表3)。

圖1 DEGs的PPI網(wǎng)絡(luò)

2.4 miRNA—核心基因調(diào)控網(wǎng)絡(luò)模塊的鑒定

本研究應(yīng)用CyTargetLinker預(yù)測(cè)可以與上述篩選出的10個(gè)核心基因相互作用的miRNA。結(jié)果顯示，在MicroCosm數(shù)據(jù)庫(kù)中有258個(gè)預(yù)測(cè)的miRNA靶點(diǎn)互作關(guān)系，在TargetScan數(shù)據(jù)庫(kù)中有1 171個(gè)預(yù)測(cè)的miRNA靶點(diǎn)互作關(guān)系，總共有875個(gè)節(jié)點(diǎn)與1 429條連線。另外，閾值(threshold)可對(duì)結(jié)果顯示的可視化調(diào)控網(wǎng)絡(luò)進(jìn)行支持?jǐn)?shù)據(jù)庫(kù)的疊加篩選，通過(guò)調(diào)設(shè)閾值可控制調(diào)控網(wǎng)絡(luò)的顯示結(jié)果，其設(shè)置范圍一般為1～3[16]。本研究將閾值設(shè)為2，結(jié)果顯示共有37個(gè)miRNA與7個(gè)靶基因存在互作關(guān)系。這些基因與預(yù)測(cè)的miRNA如表3所示。

3 討論

本研究通過(guò)對(duì)GEO數(shù)據(jù)庫(kù)中的基因芯片GSE57273進(jìn)行生物信息學(xué)分析獲得DEGs，并分析了有關(guān)這些基因富集的生物過(guò)程、細(xì)胞定位、分子功能和信號(hào)通路為機(jī)制研究提供了理論依據(jù)與研究方向，通過(guò)挖掘與核心基因互作的miRNA為PMOP的研究提供新思路。

GO與KEGG分析有助于更深入地認(rèn)識(shí)并篩選出DEGs的功能和作用。由于受到數(shù)據(jù)呈現(xiàn)空間的限制，本研究未能把這些DEGs所富集的生物學(xué)過(guò)程與信號(hào)通路全部列出。本研究由KEGG篩選出的富集通路以癌癥通路為首。從參與的DEGs來(lái)看，該通路僅包含了少部分核心基因，由于該通路涉及了大部分非核心基因，這可能導(dǎo)致它們富集的通路與PMOP不同。另外，就僅涉及的核心基因而言，現(xiàn)有關(guān)于這些基因功能的研究也存在局限性，大部分集中于腫瘤領(lǐng)域的研究，而本研究結(jié)果為這些核心基因參與PMOP提供了一定依據(jù)，有助于拓展核心基因的功能。

表3 由CyTargetLinker擴(kuò)展網(wǎng)絡(luò)分析預(yù)測(cè)的與7個(gè)核心基因互作的miRNA

本研究羅列的核心基因主要富集于RNA代謝過(guò)程、基因表達(dá)過(guò)程、PI3K-Akt信號(hào)通路和癌癥相關(guān)的信號(hào)通路等。HSP90AA1編碼的蛋白質(zhì)是一種功能類(lèi)似于同型二聚體的誘導(dǎo)型分子。HSP90AA1參與細(xì)胞的生長(zhǎng)發(fā)育過(guò)程，有研究證實(shí)使用HSP90AA1抑制劑可迅速導(dǎo)致細(xì)胞死亡，表明HSP90AA1在細(xì)胞活動(dòng)中發(fā)揮重要作用[17-18]。因而，驗(yàn)證HSP90AA1是否參與PMOP的骨偶聯(lián)的調(diào)控過(guò)程值得深入探究。EP300編碼與腺病毒E1A相關(guān)的細(xì)胞p300轉(zhuǎn)錄共激活蛋白。它與組蛋白乙酰轉(zhuǎn)移酶的功能類(lèi)似，可以通過(guò)染色質(zhì)重塑調(diào)節(jié)轉(zhuǎn)錄并且在細(xì)胞增殖和分化過(guò)程發(fā)揮重要作用。EP300在骨髓中高表達(dá)，并參與了細(xì)胞成骨分化和骨量減少的調(diào)控過(guò)程。EP300在髓核細(xì)胞中可受到骨形態(tài)發(fā)生蛋白2和骨形態(tài)發(fā)生蛋白7的調(diào)控，為EP300與PMOP的關(guān)聯(lián)提供理論支撐[19-20]。由SMARCA2基因編碼的蛋白屬于SWI/SNF家族蛋白，并且這種蛋白與果蠅的brahma蛋白高度相似。該家族蛋白具有解旋酶和ATP酶活性，其通過(guò)改變基因周?chē)娜旧|(zhì)結(jié)構(gòu)來(lái)發(fā)揮調(diào)節(jié)基因轉(zhuǎn)錄的功能。SMARCA2在卵巢、大腦等處高表達(dá)，而這類(lèi)蛋白的減少會(huì)影響間充質(zhì)干細(xì)胞的成骨分化和成脂分化平衡，進(jìn)而導(dǎo)致OP[21]。RANBP2編碼與核孔復(fù)合物免疫定位有關(guān)的RAN結(jié)合蛋白，主要在睪丸、甲狀腺、骨髓和大腦等處高表達(dá)。RAN是與核膜相關(guān)的RAS超家族的小GTP結(jié)合蛋白，它通過(guò)與蛋白質(zhì)的相互作用來(lái)調(diào)控多種細(xì)胞功能。RANBP2細(xì)胞定位為核孔復(fù)合體，并且在神經(jīng)視網(wǎng)膜中的表達(dá)非常豐富[22]。ASH1L編碼轉(zhuǎn)錄激活因子的三空腔結(jié)構(gòu)蛋白質(zhì)，它同樣在睪丸、甲狀腺、骨髓和大腦等處高表達(dá)。現(xiàn)有研究證實(shí)ASH1L可以參與細(xì)胞分化以及骨髓造血的調(diào)控[23-24]。EIF4E基因編碼的蛋白質(zhì)是真核生物翻譯起始因子4F復(fù)合物的組成部分。研究發(fā)現(xiàn)，EIF4E在病理狀態(tài)下的骨髓間充質(zhì)干細(xì)胞(bone marrow stromal cells，BMSCs)中表達(dá)水平會(huì)發(fā)生變化[25]。當(dāng)然，EIF4E作為一種原癌基因，其表達(dá)和激活與轉(zhuǎn)化和腫瘤發(fā)生密切相關(guān)。PTEN編碼的蛋白質(zhì)是磷脂酰肌醇-3，4，5-三磷酸3-磷酸酶，該基因被認(rèn)為是在大量癌癥中高頻率突變的腫瘤抑制因子。近年來(lái)有研究已經(jīng)發(fā)現(xiàn)PTEN/PI3K/AKT信號(hào)通路可能是調(diào)控BMSCs增殖和分化的途徑之一[26]。CNOT6L可能參與包括細(xì)胞增殖的各種細(xì)胞活動(dòng)，有研究證實(shí)CNOT6L可通過(guò)p53介導(dǎo)的信號(hào)通路途徑來(lái)調(diào)節(jié)細(xì)胞周期阻滯和衰老[27-28]。RPL7編碼一種由60S亞基組成的核糖體蛋白，該蛋白屬于核糖體蛋白L30P家族，主要在卵巢與骨髓高表達(dá)。KRAS編碼一種屬于小GTP酶超家族成員的蛋白。有研究證實(shí)KRAS可以與骨形態(tài)發(fā)生蛋白4靶向性結(jié)合[29]。

此外，本研究還擴(kuò)展了與核心基因相互作用的miRNA。與HSP90AA1相互調(diào)控的miRNA參與了類(lèi)固醇生物合成過(guò)程。與EP300相互作用的miRNA參與D-谷氨酰胺和D-谷氨酸代謝和細(xì)胞周期調(diào)控過(guò)程。與SMARCA2靶向調(diào)控的miRNA富集于調(diào)節(jié)干細(xì)胞多能性的信號(hào)通路、肌動(dòng)蛋白細(xì)胞骨架的調(diào)控與甲狀腺激素信號(hào)通路等。與RANBP2互作的miRNA功能富集于蛋白的內(nèi)質(zhì)網(wǎng)加工過(guò)程和雌激素信號(hào)通路等。hsa-miR-137、hsa-miR-219-5p、hsa-miR-548 d-3p可調(diào)控ASH1L，參與了細(xì)胞氮化合物代謝過(guò)程。hsa-miR-134、hsa-miR-935可調(diào)控KRAS，目前尚未發(fā)現(xiàn)它們所富集的信號(hào)通路，因此仍需進(jìn)一步研究探索。與EIF4E互作miRNA的功能富集于細(xì)胞周期調(diào)控和糖胺聚糖的生物合成等。

本研究預(yù)測(cè)的37個(gè)miRNA與核心基因?yàn)檠芯縋MOP提供了新的角度。然而，本研究所得到的結(jié)論尚需要更多的研究驗(yàn)證這些核心基因和miRNA在PMOP病理進(jìn)展中的特異性與可靠性。隨著未來(lái)針對(duì)這些核心基因及miRNA不斷地深入研究，它們或可應(yīng)用于PMOP病理的篩查，或?yàn)镻MOP藥物研發(fā)提供作用靶點(diǎn)，或有助于對(duì)PMOP合并腫瘤患者進(jìn)行基因診斷。此外，上述核心基因及miRNA或可作為遺傳致病因子的篩查點(diǎn)并為此定制個(gè)性化的防治措施。由于PMOP機(jī)制復(fù)雜，影響因素眾多，探尋并梳理具有邏輯關(guān)系和實(shí)際意義的基因及互作集簇，找出后續(xù)研究工作的切入點(diǎn)是需要進(jìn)一步思考的問(wèn)題。