周子墨 柳達 陳森相 覃森

中國醫(yī)科大學(xué)附屬盛京醫(yī)院骨科，遼寧 沈陽 110004

成骨分化(osteogenesis differentiation)是一種涉及到多基因和多步驟的復(fù)雜過程。骨折、骨質(zhì)疏松等骨病的發(fā)展和修復(fù)是成骨細胞和破骨細胞共同作用的結(jié)果。間充質(zhì)干細胞(mesenchymal stem cell，MSCs)可以分化形成成骨細胞，成骨細胞通過基質(zhì)分泌和鈣礦化，促進骨質(zhì)再生，起到重要的生物學(xué)作用[1]。

非編碼RNA(noncoding RNA，ncRNA)是一類很少甚至沒有編碼潛能的RNA。在過去的研究中，其被認為是一類轉(zhuǎn)錄過程的“噪音”。隨著高通量測序技術(shù)的發(fā)展，越來越多的noncoding RNA被識別靶向在mRNA表達的調(diào)控通路中，在調(diào)控核染色質(zhì)結(jié)構(gòu)和基因表達方面具有積極的意義[2]。長鏈非編碼RNA(long noncoding RNA，lncRNA)是一類長度在200～100 000 nt的noncoding RNA。lncRNA被分為5種類型：正義型(sense)、反義型(antisense)、雙向型(bidirectional)、內(nèi)含型(intronic)、基因間型(intergenic)，它們在功能上存在一定差異[3]。MicroRNA(miRNA)作為一類調(diào)控性RNA，可以通過多種方式介導(dǎo)基因表達，這些途徑可以被上游的lncRNA靶向調(diào)控。目前，大量研究證明lncRNA可通過內(nèi)源性競爭性RNA(ceRNA)學(xué)說靶向miRNA分子，形成miRNA海綿(Sponge)調(diào)節(jié)下游基因表達[4]。高通量測序技術(shù)可以通過基因芯片得到相應(yīng)的差異基因，被廣泛運用于基因組學(xué)和轉(zhuǎn)錄組學(xué)研究。Lim等[5]利用基因芯片，獲得大量DEGs，預(yù)測了miRNA-mRNA潛在通路。

本研究通過對數(shù)據(jù)庫中3個不同基因芯片對MSCs的成骨分化進行測序的結(jié)果進行分析，得到了差異表達的RNA及DEGs，通過靶點的互相預(yù)測，最終分析得到了差異較大的lncRNA-miRNA-mRNA互作網(wǎng)絡(luò)，同時構(gòu)建了PPI網(wǎng)絡(luò)，預(yù)測關(guān)鍵基因，并將DEGs及關(guān)鍵基因進行功能富集分析和通路富集分析，為進一步的機制研究提供理論依據(jù)。

1 材料與方法

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

從GEO公共數(shù)據(jù)庫(https://www.ncbi.nlm.nih.gov/geo/)下載基因表達芯片GSE89330、GSE72429、GSE74837。GSE89330采用樣本為成骨誘導(dǎo)分化14 d MSCs及分化的MSCs(對照組)，平臺為GPL16956，Agilent-045997 Arraystar human lncRNA microarray V3 (Probe Name Version)；GSE72429采用樣本為3個不同細胞系滑膜源性干細胞(synovial membrane-derived MSCs，SM-MSCs)和4個不同細胞系脂肪源性干細胞(adipose-derived stem cells，ADSC)來源的已分化成骨細胞，并與來自同一細胞系的未分化細胞進行比較。平臺為GPL16770，Agilent-031181 Unrestricted_Human_miRNA_V16.0_Microarray (miRBase release 16.0 miRNA ID version)。GSE74837樣本為用骨生成誘導(dǎo)方法誘導(dǎo)的2 d和5 d MSCs，對照組為未進行誘導(dǎo)的人骨髓間充質(zhì)干細胞。平臺為GPL13915 3 D-Gene Human Oligo chip 25k V2.1。

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

使用GEO2R(https://www.ncbi.nlm.nih.gov/geo/geo2r/)在線分析并獲得差異基因(differentially expressed genes，DEGs)。GEO2R是美國國立生物技術(shù)信息中心(National Center for Biotechnology Information，NCBI)的GEO數(shù)據(jù)庫自帶公共在線分析工具[6]，它將GEO數(shù)據(jù)進行l(wèi)imma差異表達分析，得出logFC、P.Value等數(shù)據(jù)。本研究將GSE74837得到的數(shù)據(jù)取P<0.05且∣logFC∣>2為DEGs；GSE72429得到的數(shù)據(jù)取P<0.05且∣logFC∣>1為DEmiRNAs；GSE89330得到的數(shù)據(jù)取P<0.05且∣logFC∣>2為DElncRNAs[7]。將得到的數(shù)據(jù)用火山圖和熱圖表達，火山圖使用R語言軟件(版本1.3.1093)R包ggplot2繪制，熱圖使用TBtools軟件(版本1.064)繪制。

1.3 基因功能富集分析

GO(gene ontology)是常用的DEGs功能作用分析方法，其可以通過注釋基因或其產(chǎn)物、識別基因芯片數(shù)據(jù)而分析得到生物學(xué)特征結(jié)果。GO富集按照生物途徑(biological process，BP)、細胞定位(cellular component，CC)、分子功能(molecular function，MF)對基因進行注釋和分類。KEGG分析(Kyoto Encyclopedia of Genes and Genomes analysis)是另一種常用的基因通路分析方法，可以將基因注釋和富集在信號通路上進行研究[8-9]。本研究應(yīng)用DAVID (Database for Annotation，Visualization and Integrated Discovery)在線數(shù)據(jù)庫分析工具進行分析，獲得所需的GO和KEGG數(shù)據(jù)。本研究使用的DAVID數(shù)據(jù)庫版本6.8，地址為https://david.ncifcrf.gov/，由美國國立變態(tài)反應(yīng)與傳染病研究所支持[10]。隨后使用Fisher Exact 或EASE Score 統(tǒng)計方法，以GO數(shù)據(jù)各項P<0.05且FDR<0.05為篩選條件，KEGG<0.05為篩選條件進行分析篩選[11-12]。

1.4 內(nèi)源性競爭性lncRNA-miRNA-mRNA互作用網(wǎng)絡(luò)構(gòu)建

內(nèi)源性競爭性RNA(competing endogenous RNA，ceRNA)是以miRNA海綿為基礎(chǔ)，lncRNA可直接與miRNAs相互作用并調(diào)控其活性的假設(shè)而構(gòu)建的。基于這一假設(shè)，通過3個步驟建立了lncRNA-miRNA-mRNA網(wǎng)絡(luò)：(1)篩選潛在可能的lncRNA、miRNA和mRNA。本研究以差異基因作為潛在可能靶點；(2)利用miRcode在線預(yù)測工具(http://www.mircode.org)根據(jù)mRNA預(yù)測潛在的miRNA靶點；(3)利用DIANA在線預(yù)測工具LncBase版本2.0(http://carolina.imis.athena-innovation.gr/diana_tools/web/index.php)進行基于miRNA的潛在lncRNA靶點預(yù)測，篩選方法為臨界值=0.7，細胞類型為hMSCs-BM，組織為骨組織，種類為干細胞[13]。最后，從三者中篩選出符合條件的RNA構(gòu)建ceRNA網(wǎng)絡(luò)。使用Cytoscape(version 3.5.1)可視化lncRNA-miRNA-mRNA網(wǎng)絡(luò)。

1.5 蛋白互作網(wǎng)絡(luò)構(gòu)建和分析

STRING (Search Tool for the Retrieval of Interacting Genes)是一款在線分析工具，可以用來構(gòu)建蛋白互作網(wǎng)絡(luò)(protein-protein interaction network，PPI network)和進行在線分析[14]。在“Creative Commons BY 4.0”中，數(shù)據(jù)可以開放獲取。將篩選的所有DEGs導(dǎo)入STRING(版本11.0，https://string-db.org/)在線分析軟件，設(shè)置置信度(confidence score)≥0.4，互作最大值(maximum number of interactors)=0，得到相應(yīng)數(shù)據(jù)[15]，并導(dǎo)入Cytoscape(版本3.5.1)進行MCODE (Molecular Complex Detection，版本1.4.2)分析，設(shè)置參數(shù)為degree cutoff=2，node score cutoff=0.2，k-core=2，max.depth=100，從PPI中篩選出連接最為緊密的集簇，并將score設(shè)置為>4，得到兩個最為緊密的集簇[16]。將得到的結(jié)果通過ClueGO+CluePedia進行GO富集分析。

1.6 關(guān)鍵基因篩選

將PPI網(wǎng)絡(luò)導(dǎo)入Cytoscape(版本3.5.1)軟件進行關(guān)鍵基因(hub gene)篩選，使用cytoHubba(版本0.1)進行分析及可視化，設(shè)置篩選條件hubba nodes ranked為degree=10，得到10個關(guān)鍵基因，并進行可視化處理[17]。

2 結(jié)果

2.1 差異表達RNA的篩選和鑒定

本研究選用基因表達芯片GSE89330、GSE72429、GSE74837，經(jīng)GEO2R初步分析和條件篩選，共獲得186個DEGs，包含81個下調(diào)基因和105個上調(diào)基因；89個DEmiRNA，包括25個下調(diào)miRNA和64個上調(diào)miRNA；441個DElncRNA，包括205個下調(diào)lncRNA和236個上調(diào)lncRNA。通過火山圖和熱圖進行可視化，得到明顯的差異表達(圖1)。

圖1 DEGs、DElncRNA、DEmiRNA的篩選和差異表達A：DEGs火山圖，紅色點代表篩選出的上調(diào)DEGs，藍色點代表篩選出的下調(diào)DEGs；B：miRNA火山圖，紅色點代表篩選出的上調(diào)DEmiRNA，藍色點代表篩選出的下調(diào)DEmiRNA；C：DElncRNA火山圖，紅色點代表篩選出的上調(diào)DElncRNA，藍色點代表篩選出的下調(diào)DElncRNA；D：GSE74837中DEGs表達熱圖；E：GSE72429中DEmiRNA表達熱圖；F：GSE89330中DElncRNA表達熱圖。Fig.1 The screen and differently expression genes in DEGs, DElncRNA, and DEmiRNAA: In the DEGs volcano map, the red dot represents the screened up-regulated DEGs, and the blue dot represents the screened down-regulated DEGs; B: MiRNA volcano map, red dot represents the screened up-regulated demirna, and blue dot represents the screened down-regulated demirna; C: In the volcano diagram of delncrna, the red dot represents the screened up-regulated delncrna, and the blue dot represents the screened down-regulated delncrna; D: Heat map of DEGs expression in gse74837; E: Heat map of demirna expression in gse72429; F: Heat map of delncrna expression in gse89330.

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

根據(jù)GO分析和KEGG分析，以P<0.05，F(xiàn)DR<0.05為篩選條件，并按照計數(shù)值從大到小排列，將BP、CC與MF類別中排名在前10的富集結(jié)果用柱狀圖表示，將KEGG分析結(jié)果用氣泡圖表示(圖2)。分析發(fā)現(xiàn)，在生物學(xué)途徑上，DEGs主要富集在炎癥反應(yīng)和血管生成和骨骼系統(tǒng)發(fā)育中，促進MSCs成骨分化。在細胞成分中，DEGs大量富集在細胞膜和胞外體以及細胞外間隙。這提示DEGs表達分泌蛋白到胞外介導(dǎo)MSCs成骨分化。在分子作用上，生長因子活化、肝素結(jié)合等與細胞內(nèi)外代謝、細胞生長和凋亡相關(guān)途徑均被富集(圖2A)。在KEGG分析中，發(fā)現(xiàn)DEGs主要富集在礦物質(zhì)吸收、TNF信號通路、malaria等通路中(圖2B)。

圖2 DEGs GO富集分析與KEGG分析A：DEGs的GO富集分析；B：DEGs的KEGG分析。Fig.2 The GO enrichment analysis and KEGG analysis of DEGsA: Go enrichment analysis of DEGs; B: KEGG analysis of DEGs

2.3 內(nèi)源性競爭性lncRNA-miRNA-mRNA互作用網(wǎng)絡(luò)構(gòu)建和分析

將篩選得到的DERNA進行比對和匹配，最終得到內(nèi)源性競爭性lncRNA-miRNA-mRNA互作用網(wǎng)

絡(luò)，并通過logFC值判斷上下調(diào)，得到可視化網(wǎng)絡(luò)(圖3)。LncRNA SNAI3-AS1、TEX41、NR2F1-AS1和下游潛在靶點miRNA hsa-miR-1207-5p連接，DANCR和下游hsa-miR-1275連接，ARHGAP22-IT1和hsa-miR-224-5p連接，MEG8、TPT1-AS1、SNHG5和hsa-miR-335-3p，SNHG1和hsa-miR-1915-3p連接，LINC01503和hsa-miR-629-5p連接，LINC00472和hsa-miR-22-3p連接，調(diào)節(jié)下游DGEs表達。

圖3 內(nèi)源性競爭性lncRNA-miRNA-mRNA互作用網(wǎng)絡(luò)Fig.3 LncRNA-miRNA-mRNA network based on ceRNAs theory注：圖中綠色方形代表lncRNA，綠色菱形代表miRNA，圓形代表DEGs，紅色代表上調(diào)DEGs，綠色代表下調(diào)DEGs，顏色越深，∣logFC∣值越大。

2.4 蛋白互作網(wǎng)絡(luò)構(gòu)建及集簇模塊篩選

通過STRING和Cytoscape軟件分析與繪制，得出PPI網(wǎng)絡(luò)(圖4A)。紅色點代表上調(diào)DEGs，綠色代表下調(diào)DEGs，顏色越深，∣logFC∣值越高。通過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，得到兩個集簇模塊，模塊1包含14個節(jié)點(node)，75條連線(edge)，score=11.538(圖4B)；模塊2包含6個節(jié)點(node)，15條連線(edge)，score=6(圖4C)。富集分析結(jié)果顯示，模塊1中DEGs主要富集在白細胞血管內(nèi)皮粘附(leukocyte adhesion to vascular endothelial cell)和中性粒細胞趨化(neutrophil chemotaxis)中(表1)。KEGG分析提示其主要富集于Malaria通路和類風(fēng)濕性關(guān)節(jié)炎(rheumatoid arthritis)相關(guān)通路中。模塊2中，DEGs主要富集在血管直徑的正向調(diào)節(jié)(positive regulation of blood vessel diameter)中，KEGG分析中，主要富集在炎癥介質(zhì)對TRP通道的調(diào)節(jié)(inflammatory mediator regulation of TRP channels)通路中。

表1 蛋白-蛋白互作網(wǎng)絡(luò)集簇DEGs的GO分析和KEGG分析Table 1 GO enrichment analysis and KEGG analysis in PPI network and significant clusters

圖4 PPI網(wǎng)絡(luò)及集簇模塊分析A：PPI網(wǎng)絡(luò)構(gòu)建，紅色圓形代表上調(diào)的表達，綠色代表下調(diào)的表達，灰色線條代表相互連接；B：集簇模塊(score=11.538)，包含14個節(jié)點(node)，75條連線(edge)；C：集簇模塊(score=6)，包含6個節(jié)點(node)，15條連線(edge)Fig.4 PPI network and significant clusters analysisA: PPI network construction, red circle represents up-regulated expression, green represents down-regulated expression, and gray line represents interconnection; B: Cluster module (score=11.538), including 14 nodes and 75 edges; C: Cluster module (score=6), including 6 nodes and 15 edges.

2.5 關(guān)鍵基因篩選

使用cytoHubba(版本0.1)進行篩選，獲得了10個關(guān)鍵基因(hub gene)(圖5A)，顏色越深，代表ranked score排名越靠前。10個hub gene分別為IL6、CXCL12、CXCL8、CCL2、HGF、LEP、VCAM1、CXCL1、SAA1、FOS。對這10個hub gene進行GO分析，發(fā)現(xiàn)其主要富集在趨化因子活化、白細胞血管粘附和單核細胞趨化中(圖5B、5D)。KEGG分析顯示hub gene主要以Malaria相關(guān)通路和白血球粘附血管內(nèi)皮相關(guān)通路為主(圖5C)。

圖5 基于degree方法在PPI網(wǎng)絡(luò)中的前10個關(guān)鍵基因A：關(guān)鍵基因篩選及互作網(wǎng)絡(luò)，顏色越深，代表ranked score排名越靠前；B：關(guān)鍵基因GO富集分析及法分類；C：KEGG富集分析及分類；D：關(guān)鍵基因GO富集直條圖，縱坐標(biāo)代表GO通路，橫坐標(biāo)代表富集的基因數(shù)目。Fig.5 Top 10 hub genes in PPI network ranked by Degree methodA: Key gene screening and interaction network, the darker the color, the higher the ranking of ranked score; B: Go enrichment analysis and classification of key genes; C: KEGG enrichment analysis and classification; D: Go enrichment bar graph of key genes, the ordinate represents go pathway, and the abscissa represents the number of enriched genes.

3 討論

lncRNA和miRNA可以調(diào)節(jié)成骨分化中多種信號通路增強或抑制成骨分化能力[18-19]。本研究對GSE89330、GSE72429、GSE74837三組數(shù)據(jù)的分析和篩選，得到了lncRNA-miRNA-mRNA互作網(wǎng)絡(luò)。Fang等[20]發(fā)現(xiàn)，DANCR可以通過調(diào)節(jié)滑膜液源性MSCs的miR-1275/MMP-13軸，誘導(dǎo)軟骨形成。Chao等[21]通過對激素性股骨頭壞死的患者miR-1207-5p的分析發(fā)現(xiàn)，患者miR-1207-5p水平顯著升高，可能靶向抑制VEGF表達，抑制骨再生。本實驗結(jié)果中同樣獲得了DANCR、miR-1207-5p、miR-1275等RNA的差異表達。這些潛在靶點可以為分子層面研究提供更多的理論依據(jù)。

通過對DGEs的GO分析及KEGG分析，發(fā)現(xiàn)DEGs富集在炎癥反應(yīng)、血管生成和骨骼系統(tǒng)發(fā)育中。在KEGG分析中，發(fā)現(xiàn)DEGs主要富集在礦物質(zhì)吸收、TNF信號通路、malaria等通路中。在構(gòu)建PPI互作網(wǎng)絡(luò)后，通過集簇篩選和hub gene篩選得到最密集集簇和關(guān)鍵基因。通過再次GO分析和KEGG分析，將富集的范圍縮小，得到最終的富集通路結(jié)果。關(guān)鍵基因主要富集在白細胞系帶或滾動(GO:0050901)、白細胞血管內(nèi)皮粘附(GO:0061756)及趨化因子活化(GO:0008009)等免疫相關(guān)通路中。成骨細胞可以產(chǎn)生免疫相關(guān)性蛋白和趨化因子，誘導(dǎo)B細胞的成熟。趨化因子CXCL12被證明在此過程中被上調(diào)，體內(nèi)實驗已經(jīng)證明給予CXCL12后可以促進MSCs的聚集[22-24]。在本實驗中，篩選出的關(guān)鍵基因中，CXCL12被證明顯著下調(diào)，證明大量的趨化因子CXCL12靶向作用在MSCs細胞上。

綜上，本研究通過構(gòu)建lncRNA-miRNA-mRNA的互作網(wǎng)絡(luò)，顯示了SNHG1/miR-1915-3p/CXCL12、LINC00472/miR-22-3p/HSD11B1、LINC01503/miR-629-3p/ZBTB16、DANCR/miR-1275/MMP7等具有潛在靶向的通路。分析獲得的相關(guān)差異基因IL6、CXCL12、CXCL8、CCL2、HGF、LEP、VCAM1、CXCL1、SAA1、FOS和ncRNA均在成骨分化中具有較高的組織特異性和功能特異性。