成骨分化內源性競爭性lncRNA-miRNA-mRNA網絡與核心基因的研究

周子墨 柳達 陳森相 覃森

中國醫科大學附屬盛京醫院骨科，遼寧 沈陽 110004

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

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

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

1 材料與方法

1.1 數據收集

從GEO公共數據庫(https://www.ncbi.nlm.nih.gov/geo/)下載基因表達芯片GSE89330、GSE72429、GSE74837。GSE89330采用樣本為成骨誘導分化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樣本為用骨生成誘導方法誘導的2 d和5 d MSCs，對照組為未進行誘導的人骨髓間充質干細胞。平臺為GPL13915 3 D-Gene Human Oligo chip 25k V2.1。

1.2 差異基因的數據分析

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

1.3 基因功能富集分析

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

1.4 內源性競爭性lncRNA-miRNA-mRNA互作用網絡構建

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

1.5 蛋白互作網絡構建和分析

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

1.6 關鍵基因篩選

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

2 結果

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

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

圖1 DEGs、DElncRNA、DEmiRNA的篩選和差異表達A：DEGs火山圖，紅色點代表篩選出的上調DEGs，藍色點代表篩選出的下調DEGs；B：miRNA火山圖，紅色點代表篩選出的上調DEmiRNA，藍色點代表篩選出的下調DEmiRNA；C：DElncRNA火山圖，紅色點代表篩選出的上調DElncRNA，藍色點代表篩選出的下調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信號通路富集分析

根據GO分析和KEGG分析，以P<0.05，FDR<0.05為篩選條件，并按照計數值從大到小排列，將BP、CC與MF類別中排名在前10的富集結果用柱狀圖表示，將KEGG分析結果用氣泡圖表示(圖2)。分析發現，在生物學途徑上，DEGs主要富集在炎癥反應和血管生成和骨骼系統發育中，促進MSCs成骨分化。在細胞成分中，DEGs大量富集在細胞膜和胞外體以及細胞外間隙。這提示DEGs表達分泌蛋白到胞外介導MSCs成骨分化。在分子作用上，生長因子活化、肝素結合等與細胞內外代謝、細胞生長和凋亡相關途徑均被富集(圖2A)。在KEGG分析中，發現DEGs主要富集在礦物質吸收、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 內源性競爭性lncRNA-miRNA-mRNA互作用網絡構建和分析

將篩選得到的DERNA進行比對和匹配，最終得到內源性競爭性lncRNA-miRNA-mRNA互作用網

絡，并通過logFC值判斷上下調，得到可視化網絡(圖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連接，調節下游DGEs表達。

圖3 內源性競爭性lncRNA-miRNA-mRNA互作用網絡Fig.3 LncRNA-miRNA-mRNA network based on ceRNAs theory注：圖中綠色方形代表lncRNA，綠色菱形代表miRNA，圓形代表DEGs，紅色代表上調DEGs，綠色代表下調DEGs，顏色越深，∣logFC∣值越大。

2.4 蛋白互作網絡構建及集簇模塊篩選

通過STRING和Cytoscape軟件分析與繪制，得出PPI網絡(圖4A)。紅色點代表上調DEGs，綠色代表下調DEGs，顏色越深，∣logFC∣值越高。通過MCODE(Molecular Complex Detection，版本1.4.2)，設置參數為degree cutoff=2，node score cutoff=0.2，k-core=2，max.depth=100，score>4，得到兩個集簇模塊，模塊1包含14個節點(node)，75條連線(edge)，score=11.538(圖4B)；模塊2包含6個節點(node)，15條連線(edge)，score=6(圖4C)。富集分析結果顯示，模塊1中DEGs主要富集在白細胞血管內皮粘附(leukocyte adhesion to vascular endothelial cell)和中性粒細胞趨化(neutrophil chemotaxis)中(表1)。KEGG分析提示其主要富集于Malaria通路和類風濕性關節炎(rheumatoid arthritis)相關通路中。模塊2中，DEGs主要富集在血管直徑的正向調節(positive regulation of blood vessel diameter)中，KEGG分析中，主要富集在炎癥介質對TRP通道的調節(inflammatory mediator regulation of TRP channels)通路中。

表1 蛋白-蛋白互作網絡集簇DEGs的GO分析和KEGG分析Table 1 GO enrichment analysis and KEGG analysis in PPI network and significant clusters

圖4 PPI網絡及集簇模塊分析A：PPI網絡構建，紅色圓形代表上調的表達，綠色代表下調的表達，灰色線條代表相互連接；B：集簇模塊(score=11.538)，包含14個節點(node)，75條連線(edge)；C：集簇模塊(score=6)，包含6個節點(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 關鍵基因篩選

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

圖5 基于degree方法在PPI網絡中的前10個關鍵基因A：關鍵基因篩選及互作網絡，顏色越深，代表ranked score排名越靠前；B：關鍵基因GO富集分析及法分類；C：KEGG富集分析及分類；D：關鍵基因GO富集直條圖，縱坐標代表GO通路，橫坐標代表富集的基因數目。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可以調節成骨分化中多種信號通路增強或抑制成骨分化能力[18-19]。本研究對GSE89330、GSE72429、GSE74837三組數據的分析和篩選，得到了lncRNA-miRNA-mRNA互作網絡。Fang等[20]發現，DANCR可以通過調節滑膜液源性MSCs的miR-1275/MMP-13軸，誘導軟骨形成。Chao等[21]通過對激素性股骨頭壞死的患者miR-1207-5p的分析發現，患者miR-1207-5p水平顯著升高，可能靶向抑制VEGF表達，抑制骨再生。本實驗結果中同樣獲得了DANCR、miR-1207-5p、miR-1275等RNA的差異表達。這些潛在靶點可以為分子層面研究提供更多的理論依據。

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

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