用生物信息學分析預測絕經后骨質疏松癥核心基因與互作miRNA的研究

柴毅 譚峰 樊巧玲

南京中醫藥大學，江蘇 南京 210046

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

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

本研究通過對基因芯片GSE57273進行生物信息學分析篩選核心基因，并預測與這些核心基因相互作用的miRNA，為PMOP建立新的科學假說以及后續更深入的研究提供依據，并為PMOP診斷以及治療藥物的研發提供較為可靠的路徑和作用靶點。

1 材料與方法

1.1 數據收集

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

1.2 差異基因的數據分析

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

1.3 基因功能富集分析

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

1.4 蛋白互作網絡與模型分析

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

1.5 miRNA-基因調控網絡的預測

CyTargetLinker可以擴展生物調控互作網絡(regulatory interaction networks，RegINs)，由荷蘭系統生物學聯合會提供支持。它涵蓋了miRNA—靶點、轉錄因子—靶點和藥物—靶點之間的互作關系。本研究下載了人類物種基因數據集(https://projects.bigcat.unimaas.nl/cytargetlinker/regins/)。選用該數據集中基于實驗驗證的miRTarBase 4.4數據庫(含20 942個RegINs)，基于預測功能的TargetScan 6.2數據庫(含511 040個RegINs)和MicroCosm 5數據庫(含541 039個RegINs)預測核心基因與miRNA之間的調控關系。

2 結果

2.1 DEGs的鑒定

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

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

根據GO的分析結果，本研究以P<0.05，FDR<0.05為篩選條件，并按照計數值從大到小排列，在BP、CC與MF類別中各選取前3項列為表1。可以看出，在生物學過程中，這些DEGs主要參與了基因表達，細胞大分子生物合成和RNA代謝過程；在細胞定位中，這些DEGs富集于核漿、細胞質基質以及粘附連接；從分子功能上看，這些DEGs具有使有機環狀化合物結合、雜環化合物結合和核酸結合的作用。根據KEGG分析結果，本研究以P<0.05為篩選條件，按照計數值從大到小將DEGs富集的信號通路列為表2。結果顯示，這些DEGs主要富集于癌癥信號通路、病毒致癌通路、粘附斑、rap1信號通路和內質網蛋白加工通路。

表1 與PMOP相關的DEGs的GO富集分析

表2 與PMOP相關的DEGs的KEGG富集分析

圖2 蛋白質互作網絡的前3個集簇模塊

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

通過STRING的PPI構建，經由Cytoscape對網絡的計算工具得出所有DEGs的連接度(degree)(圖1)。degree值表示網絡中某一基因與周圍基因的關系數量，因此degree越大代表與它相互作用關系的基因數量就越多。本研究按degree從高至低進行排序，以排名前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分析發現17組集簇，共包含523個節點(node)與2 026條連線(edge)。本研究以score為依據展示前3組集簇并分析各個集簇所富集的通路(圖2)。其中集簇模塊A富集核糖體(hsa03010：Ribosome)與mRNA監視通路(hsa03015：mRNA surveillance pathway)(P<0.05)，集簇模塊B富集泛素介導的蛋白質水解4通路(hsa04120：Ubiquitin mediated proteolysis 4)(P<0.05)，而集簇模塊C未鑒定出具有統計學意義的通路(P>0.05)(表3)。

圖1 DEGs的PPI網絡

2.4 miRNA—核心基因調控網絡模塊的鑒定

本研究應用CyTargetLinker預測可以與上述篩選出的10個核心基因相互作用的miRNA。結果顯示，在MicroCosm數據庫中有258個預測的miRNA靶點互作關系，在TargetScan數據庫中有1 171個預測的miRNA靶點互作關系，總共有875個節點與1 429條連線。另外，閾值(threshold)可對結果顯示的可視化調控網絡進行支持數據庫的疊加篩選，通過調設閾值可控制調控網絡的顯示結果，其設置范圍一般為1～3[16]。本研究將閾值設為2，結果顯示共有37個miRNA與7個靶基因存在互作關系。這些基因與預測的miRNA如表3所示。

3 討論

本研究通過對GEO數據庫中的基因芯片GSE57273進行生物信息學分析獲得DEGs，并分析了有關這些基因富集的生物過程、細胞定位、分子功能和信號通路為機制研究提供了理論依據與研究方向，通過挖掘與核心基因互作的miRNA為PMOP的研究提供新思路。

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

表3 由CyTargetLinker擴展網絡分析預測的與7個核心基因互作的miRNA

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

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

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