基于生物信息學(xué)分析老年骨質(zhì)疏松差異表達(dá)基因及m6A相關(guān)蛋白的研究

周子墨，陳佳慧，陳森相，覃森，黃瑛，柳達(dá)

（中國(guó)醫(yī)科大學(xué)附屬盛京醫(yī)院 1.骨科；2.超聲科，沈陽(yáng) 110004）

骨質(zhì)疏松是一類(lèi)常見(jiàn)的代謝性疾病，在絕經(jīng)后婦女及老年人中常見(jiàn)。骨礦物質(zhì)的流失導(dǎo)致骨密度量下降及骨微結(jié)構(gòu)失衡和退化，骨骼的脆性和骨折風(fēng)險(xiǎn)大大增加［1］。人骨髓源性間充質(zhì)干細(xì)胞（human bone marrow derived mesenchymal stem cells，hBMSCs）作為一類(lèi)多能干細(xì)胞，其成骨和成脂分化的方向選擇對(duì)骨質(zhì)疏松癥進(jìn)展有較大影響［2］。

近年來(lái)，隨著測(cè)序技術(shù)的發(fā)展，多種基因及其產(chǎn)物可以通過(guò)基因表達(dá)譜的分析和可視化與疾病相關(guān)聯(lián)［3］。在不同疾病或同一疾病不同亞型中，基因表達(dá)都有較大差異。m6A甲基化相關(guān)蛋白在多種疾病中存在一定的調(diào)控能力，在hBMSCs中也發(fā)揮著重要的作用［4］。

本研究利用生物信息學(xué)技術(shù)，分析老年骨質(zhì)疏松癥患者h(yuǎn)BMSCs的基因芯片，獲得差異表達(dá)基因（differentially expressed genes，DEGs）譜，并構(gòu)建m6A相關(guān)蛋白差異表達(dá)譜，為老年骨質(zhì)疏松癥的診治和預(yù)防提供新的思路與方向。

1 材料與方法

1.1 數(shù)據(jù)收集與篩選

從NCBI 基因表達(dá)綜合（gene expression omnibus，GEO）數(shù)據(jù)庫(kù)（http：//www.ncbi.nlm.nih.gov/geo/）下載基因表達(dá)芯片GSE35956。GSE35956采用樣本為從老年骨質(zhì)疏松患者和老年非骨質(zhì)疏松患者骨髓提取的hBMSCs，平臺(tái)為GPL570，Affymetrix Human Genome U133 Plus 2.0 Array。

1.2 DEGs數(shù)據(jù)分析

使用R軟件的Limma軟件包（版本：3.40.2）分析數(shù)據(jù)集中mRNA的差異表達(dá)情況［5］。數(shù)據(jù)選取adjustedP＜ 0.05且|logFC|＞2（logFC被定義為閾值mRNA差異表達(dá)的篩選）為DEGs。DEGs熱圖采用R軟件包pheatmap繪制。

1.3 基因功能富集分析

為進(jìn)一步篩選并確認(rèn)DEGs的潛在功能，通過(guò)基因本體論（gene ontology，GO）對(duì)數(shù)據(jù)進(jìn)行了功能富集分析。GO是一種廣泛使用的工具，通過(guò)分子功能（molecular function，MF）、生物途徑（biological pathway，BP）和細(xì)胞成分（cellular constituent，CC）來(lái)注釋具有潛在功能的基因［6］。京都基因和基因組數(shù)據(jù)庫(kù)（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析是一種用于注釋基因功能以及相關(guān)信號(hào)通路的分析方法。本研究使用R語(yǔ)言軟件中的Cluster-Profiler程序包分析DEGs的GO功能富集以及KEGG通路，氣泡圖使用R軟件包ggord進(jìn)行繪制。

1.4 m6A相關(guān)蛋白表達(dá)譜分析

候選的20種m6A相關(guān)蛋白基因來(lái)源于多種系及多種疾病的m6A相關(guān)分子及機(jī)制研究。m6A相關(guān)蛋白基因差異分析采用wilcox檢驗(yàn)，箱線圖通過(guò)R軟件包 ggplot2 進(jìn)行實(shí)現(xiàn)；熱圖通過(guò)R軟件包pheatmap進(jìn)行繪制。

1.5 蛋白質(zhì)-蛋白質(zhì)相互作用（protein-protein interaction，PPI）網(wǎng)絡(luò)構(gòu)建、集簇分析和關(guān)鍵（Hub）基因獲取

STRING（Search Tool for the Retrieval of Interacting Genes）在線分析工具可以實(shí)現(xiàn)PPI網(wǎng)絡(luò)可視化分析［7］。在“Creative Commons BY 4.0”中，數(shù)據(jù)可以開(kāi)放獲取。導(dǎo)入已篩選DEGs進(jìn)入 STRING（版本11.0，https：//string-db.org/），設(shè)置置信度＞0.4，互作最大值=0，進(jìn)行數(shù)據(jù)分析。獲得所有結(jié)果由Cytoscape（版本3.5.1）進(jìn)行可視化分析。MCODE（Molecular Complex Detection，版本1.4.2）為Cytoscape內(nèi)置分析軟件，可以從PPI中篩選出連接最為緊密的集簇。設(shè)置參數(shù)degree cutoff=2，node score cutoff=0.2，k-core=2，max.depth=100，取score值＞4，得到最為緊密的集簇。使用Cytoscape內(nèi)置cytoHubba（版本0.1）軟件進(jìn)行Hub基因分析及可視化，篩選條件hubba nodes ranked取degree=10，得到10個(gè)Hub基因。

1.6 統(tǒng)計(jì)學(xué)分析

采用R語(yǔ)言進(jìn)行數(shù)據(jù)分析。差異分析采用Wilcoxon檢驗(yàn)，P＜ 0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 DEGs的篩選和鑒定

本研究選用GEO數(shù)據(jù)集GSE35956，經(jīng)Limma分析和條件篩選，共獲得292個(gè)DEGs，包括 93個(gè)低表達(dá)基因和199個(gè)高表達(dá)基因。通過(guò)火山圖和熱圖進(jìn)行可視化，得到明顯的差異表達(dá)，DEGs個(gè)數(shù)較多，此處分別展示差異改變最大的50個(gè)上調(diào)基因和50個(gè)下調(diào)基因，見(jiàn)圖1。

圖1 DEGs熱圖及火山圖Fig.1 Heat maps and volcanic maps of differentially expressed genes

2.2 基因功能及通路富集分析

通過(guò)對(duì)DEGs的GO及KEGG分析，針對(duì)292個(gè)基因進(jìn)行了注釋?zhuān)謩e得到了上調(diào)和下調(diào)DEGs的GO富集結(jié)果及KEGG分析結(jié)果。在GO分析中，上調(diào)的DEGs主要富集在成骨過(guò)程、結(jié)締組織發(fā)展、細(xì)胞趨化、胞內(nèi)運(yùn)輸?shù)恼{(diào)節(jié)、細(xì)胞外結(jié)構(gòu)組織等；下調(diào)的DEGs主要富集在tRNA相關(guān)代謝過(guò)程、突 觸的組織等。上調(diào)和下調(diào)的DEGs均可富集在成骨過(guò)程中。KEGG通路中，高表達(dá)基因富集在PI3K-Akt信號(hào)通路以及黏著斑相關(guān)通路中，低表達(dá)基因富集在Wnt信號(hào)通路、Hippo信號(hào)通路以及一些免疫調(diào)控通路中。見(jiàn)圖2。

圖2 DEGs的GO分析及KEGG分析氣泡圖Fig.2 Bubble graphs of gene ontology and Kyoto Encyclopedia of Genes and Genomes analyses

2.3 m6A相關(guān)蛋白表達(dá)譜分析

20個(gè)候選的m6A相關(guān)蛋白中，通過(guò)差異分析得到ALKBH5、FTO、IGF2BP2、RBMX、VIRMA、YTH-DC2、YTHDF2、ZC3H13等8個(gè)基因（圖3A），并獲得基因差異分析的箱線圖（圖3B）。

圖3 m6A相關(guān)蛋白基因差異表達(dá)熱圖及差異分析Fig.3 Heat map and DEGs analysis of m6A related protein

2.4 PPI網(wǎng)絡(luò)及集簇分析和Hub基因篩選

將STRING分析的結(jié)果導(dǎo)入Cytoscape軟件進(jìn)行進(jìn)一步分析與繪制，得到PPI網(wǎng)絡(luò)。通過(guò)MCODE（Molecular Complex Detection，版本1.4.2），設(shè)置參數(shù)為degree cutoff=2，node score cutoff=0.2，k-core=2，max.depth=100，score＞4，得到2個(gè)集簇模塊，模塊1包含13個(gè)節(jié)點(diǎn)（node），37條連線（edge），score=6.167（圖4A）；模塊2包含5個(gè)節(jié)點(diǎn)（node），10條連線（edge），score=5（圖4B）。cytoHubba（版本0.1）進(jìn)行Hub基因篩選，獲得了10個(gè)Hub基因（圖4C），顏色越深，代表ranked score排名越靠前。10個(gè)Hub基因分別為CCR5、TAC1、GNGT1、RAC2、VAMP2、PXN、CSF1、TRIP10、RAB11B、ADRA2A。

圖4 PPI網(wǎng)絡(luò)、集簇和Hub基因Fig.4 Significant clusters analysis，and Hub gene

3 討論

近年來(lái)，隨著研究的不斷深入，越來(lái)越多的證據(jù)表明老年骨質(zhì)疏松的進(jìn)展與成骨分化和破骨作用失衡有關(guān)。這一過(guò)程中，hBMSCs發(fā)揮著重要的作用。本研究通過(guò)生物信息學(xué)技術(shù)，研究了老年骨質(zhì)疏松患者和非骨質(zhì)疏松患者的hBMSCs DEGs譜，通過(guò)表達(dá)譜的GO富集分析及KEGG通路分析，發(fā)現(xiàn)上調(diào)及下調(diào)的差異基因均可以富集在骨形成相關(guān)通路及細(xì)胞生長(zhǎng)分化相關(guān)通路中；而骨小梁形成過(guò)程富集在下調(diào)的差異基因中。骨小梁的形成是骨皮質(zhì)在松質(zhì)骨內(nèi)的延伸部分，對(duì)于成骨細(xì)胞及血管形成具有重要作用［8］。骨小梁的減少會(huì)極大地影響成骨細(xì)胞分化和骨再生，導(dǎo)致骨的脆性增加［9］。GO富集的結(jié)果顯示在老年骨質(zhì)疏松中，hBMSCs的差異基因?qū)巧删哂兄匾恼{(diào)控作用。KEGG通路結(jié)果顯示，PI3K-Akt信號(hào)通路和Wnt通路等被高度富集。PI3K-Akt信號(hào)通路可以通過(guò)調(diào)控下游不同的分子，促進(jìn)間充質(zhì)干細(xì)胞的成骨分化［10-12］，但TAN等［13］研究發(fā)現(xiàn)，TAZ通過(guò)抑制PI3K-Akt信號(hào)通路促進(jìn)骨質(zhì)疏松大鼠模型間充質(zhì)干細(xì)胞成骨分化。這提示PI3KAkt信號(hào)通路的復(fù)雜調(diào)控模式具有雙向性。此外，Wnt通路及Hippo信號(hào)通路已被證明可以促進(jìn)成骨分化［14-15］，其相關(guān)基因在本研究中低表達(dá)，表明成骨分化過(guò)程被明顯抑制。

m6A作為真核生物RNA中最豐富的甲基化修飾，在間充質(zhì)干細(xì)胞成骨分化的表觀遺傳上具有重要的作用［6，16］。GUO等［17］研究發(fā)現(xiàn)，F(xiàn)TO作為去甲基化酶，不僅可以對(duì)脂肪調(diào)控，還可以在骨質(zhì)疏松中發(fā)揮調(diào)控作用。本研究通過(guò)對(duì)m6A甲基化蛋白的篩選和差異分析，獲得了ALKBH5、FTO、IGF2BP2、RBMX、VIRMA、YTHDC2、YTHDF2、ZC3H13等8個(gè)基因的表達(dá)差異，包括了m6A甲基化蛋白的各種類(lèi)型。這為進(jìn)一步研究老年骨質(zhì)疏松m6A甲基化修飾提供了思路。

通過(guò)PPI網(wǎng)絡(luò)的構(gòu)建、集簇分析和Hub基因篩選得到最密集集簇和Hub基因（CCR5、TAC1、GNGT1、RAC2、VAMP2、PXN、CSF1、TRIP10、RAB11B、ADRA2A）。這些Hub基因中，CCR5、TAC1、VAMP2、CSF1等均和hBMSCs的免疫及炎癥反應(yīng)相關(guān)［18］。而骨質(zhì)疏松的發(fā)生發(fā)展過(guò)程中，炎癥反應(yīng)及其趨化因子對(duì)間充質(zhì)干細(xì)胞的遷移和歸巢起到趨化作用［19］。這也提示了骨質(zhì)疏松和免疫系統(tǒng)的相關(guān)性。Hub基因的獲得為研究間充質(zhì)干細(xì)胞在老年骨質(zhì)疏松中蛋白互相作用調(diào)控以及信號(hào)通路提供了條件。

綜上所述，本研究構(gòu)建了老年骨質(zhì)疏松hBMSCs DEGs譜，并分析了具有潛在靶向的信號(hào)通路及細(xì)胞生物學(xué)過(guò)程，還得到了m6A相關(guān)mRNA表達(dá)的差異性特征。同時(shí)，分析獲得的相關(guān)差異基因在老年骨質(zhì)疏松的hBMSCs成骨分化中具有較高的組織特異性和功能特異性。