肖賀方 連笑宇 車凌儀 強(qiáng)碩 鄭稼*

1.鄭州大學(xué)人民醫(yī)院 河南省人民醫(yī)院骨科，河南 鄭州 450003 2.鄭州大學(xué)第一附屬醫(yī)院婦產(chǎn)科，河南 鄭州 450003

絕經(jīng)后骨質(zhì)疏松癥(postmenopausal osteoporosis，PMOP)，是一種無癥狀全身漸進(jìn)性疾病，特點(diǎn)是骨密度低、骨脆性增加，從而骨折的發(fā)生率增高。疾病的關(guān)鍵環(huán)節(jié)是因絕經(jīng)期卵泡刺激素的量增加和雌激素缺乏，從而導(dǎo)致嚴(yán)重的骨量丟失[1]，研究表明成骨細(xì)胞和破骨細(xì)胞的功能對(duì)骨骼的修復(fù)和重塑至關(guān)重要[2]。目前，常用的治療方法大多使用利塞膦酸鹽，唑來膦酸等；最近新增的一些生物療法如選擇性雌激素受體拮抗劑[3]。使用雙膦酸化合物治療會(huì)產(chǎn)生許多不良反應(yīng)，包括肌肉骨骼疼痛、下頜骨折等問題[4]。在美國，每年因骨質(zhì)疏松癥，可導(dǎo)致約150萬例骨折的發(fā)生[5]。然而，作為一個(gè)全球性的問題，目前仍然缺乏有效的診斷和治療。因此，研究PMOP形成的分子機(jī)制具有十分重要的作用。在之前的文章中，一些與PMOP形成相關(guān)的潛在的差異表達(dá)基因、中心基因和信號(hào)通路尚未得到確認(rèn)，此外，沒有利用GEO數(shù)據(jù)庫，來研究此問題。

1 材料和方法

1.1 數(shù)據(jù)來源

GEO(http：//www.ncbi.nlm.nih.gov/geo)[6]是一個(gè)公共基因組學(xué)數(shù)據(jù)存儲(chǔ)庫。數(shù)據(jù)集GSE56116下載于 GPL4133 平臺(tái)，數(shù)據(jù)集有13個(gè)樣本，包含3個(gè)健康的和10個(gè)疾病的樣本。

1.2 差異分析

使用了R軟件(版本3.6.1；https：//www.r_project.org/)和對(duì)應(yīng)的R包進(jìn)行分析，在 R 軟件中執(zhí)行以下參數(shù)，包括分位數(shù)、RMA、Median polish和 pmonly。再用LIMMA包進(jìn)行差異分析，篩選標(biāo)準(zhǔn)：log2 fold change(FC)>1 和Pvalue<0.05。

1.3 富集分析

通過數(shù)據(jù)庫DAVID用于探索分析差異基因的功能[7]。GO(gene ontology)數(shù)據(jù)庫里面包括有描述基因基本特征的詞匯或結(jié)構(gòu)[8]。KEGG(Kyoto Encyclopedia of Genes and Genomes)數(shù)據(jù)庫從基因組、系統(tǒng)性和化學(xué)功能通路等方面綜合信息[9]。使用DAVID軟件，把所有的差異基因進(jìn)行功能富集分析。

1.4 構(gòu)建蛋白質(zhì)互作網(wǎng)絡(luò)(PPI)和篩選關(guān)鍵基因

PPI網(wǎng)絡(luò)通過尋找互作基因來進(jìn)行預(yù)測(STRING；http://string-db.org)(版本12.0)[10]。PPI網(wǎng)絡(luò)由Cytoscape繪制，最重要的模塊由MCODE插件進(jìn)行確認(rèn)。標(biāo)準(zhǔn)：degree cut-off=2,MCODE scores >4,Max depth=100,node score cut-off=0.3,和 k-score=2。然后，使用Hubba插件挑選中心基因。使用 DAVID對(duì)網(wǎng)絡(luò)中關(guān)鍵模塊中的基因進(jìn)行KEGG和GO分析。使用五種方法(degree、BottleNeck、Radiality centrality、Stress centrality、Closeness centrality)來排列和篩選中心基因[11]。

1.5 用CIBERSORT 算法分析免疫細(xì)胞的浸潤

利用CIBERSORT算法對(duì)獲得的數(shù)據(jù)進(jìn)行了分析，并得到了22種免疫細(xì)胞的比例。基于P值 <0.05 篩選，并計(jì)算樣本中每個(gè)免疫細(xì)胞的百分比。通過R軟件中的vioplot包，分析兩組免疫細(xì)胞的不同免疫滲透水平。

1.6 在數(shù)據(jù)庫中進(jìn)行驗(yàn)證

GSE7429 數(shù)據(jù)集是從GEO數(shù)據(jù)庫獲得，使用此數(shù)據(jù)集驗(yàn)證得到的關(guān)鍵基因。

2 結(jié)果

2.1 確定差異基因

從數(shù)據(jù)中篩選出了308個(gè)差異基因，其中110個(gè)下調(diào)，198個(gè)上調(diào)。圖1和圖2分別顯示了所有DEG的火山圖和前50個(gè)下調(diào)基因和50個(gè)上調(diào)基因的熱圖。

注：紅色的差異基因代表log2FC>1 并且P<0.05，藍(lán)色的差異基因代表log2FC<-1 并且 P<0.05。圖1 所有差異基因的火山圖Fig.1 The volcano plot of all DEGs

圖2 紅色表示熱圖中前50上調(diào)的差異基因，藍(lán)色表示前50個(gè)下調(diào)的差異基因Fig.2 The heat map of the top 50 upgrade genes in red and 50 downgrade genes in blue

2.2 GO和KEGG的富集分析

使用 DAVID 對(duì)差異基因進(jìn)行 GO 和KEGG 富集分析，如圖 3 所示，分別顯示了在 GO 分析中的前10個(gè)富集的部分。生物過程(BP)分析結(jié)果表明，DEG在先天免疫反應(yīng)(innate immune response)、炎癥反應(yīng)(inflammatory response)中明顯富集(見圖3a)；分子功能 (MF)分析表明，DEG 在碳水化合物結(jié)合(carbohydrate binding)、受體活性(receptor activity)、轉(zhuǎn)錄活化活性(transcriptional activator activity)、受體結(jié)合 (receptor binding)富集顯著(見圖3b)；對(duì)于細(xì)胞組分(CC)分析表明，DEG 主要富集于等離子膜(plasma membrane)，等離子膜的整體成分(integral component of plasma membrane)，見圖3c。

注：a 生物過程；b 分子功能；c 細(xì)胞組分；d 上、下調(diào)的差異基因的KEGG通路分析。圖3 差異基因的GO和KEGG的富集分析Fig.3 The Gene Ontology (GO) analysis of differentially expressed genes

2.3 KEGG通路分析

KEGG通路富集分析表明，DEG主要在破骨細(xì)胞分化(osteoclast differentiation)、造血細(xì)胞系(hematopoietic cell lineage)、細(xì)胞粘附分子(cell adhesion molecules CAMs)、甲狀腺素合成(thyroid hormone synthesis)途徑中富集顯著。在這些通路中破骨細(xì)胞分化的通路占主要部分。KEGG通路富集分析結(jié)果見圖3 d。

2.4 構(gòu)建PPI網(wǎng)絡(luò)

節(jié)點(diǎn)表示 DEG 和邊緣表示 DEG之間的交互連接，根據(jù)STRING數(shù)據(jù)庫中的信息，我們使用Cytoscape中的插件MCODE在PPI網(wǎng)絡(luò)中查找165個(gè)節(jié)點(diǎn)和339個(gè)邊(見圖4a)和得到了最顯著的模塊(見圖4b)。

圖4 a 使用Cytoscape構(gòu)建的差異基因的蛋白質(zhì)互作網(wǎng)絡(luò)；b 從蛋白質(zhì)互作網(wǎng)絡(luò)選取的最重要的模塊Fig.4 (a) The PPI network of DEGs was constructed using Cytoscape.(b) The most significant module was obtained from PPI network

使用DAVID對(duì)最顯著模塊中基因進(jìn)行了富集分析。結(jié)果表明，該模塊中的基因主要在GO的炎癥反應(yīng)(inflammatory response)、先天免疫反應(yīng)(innate immune response)和血漿膜(plasma membrane)中富集顯著，而在KEGG中主要在Fc伽馬R介導(dǎo)的噬菌體(Fc gamma R-mediated phagocytosis)、骨細(xì)胞分化(osteoclast differentiation)和化學(xué)因子信號(hào)(chemokine signaling)等通路中富集，見表1。

表1 對(duì)關(guān)鍵模塊的差異基因進(jìn)行GO和KEGG富集分析的前10項(xiàng)富集結(jié)果Table 1 The top 10 Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched for the genes involved in significant module

2.5 確定關(guān)鍵基因

使用Cytoscape中的cytoHubba插件來確定關(guān)鍵基因。根據(jù)cytoHubba的5種分類方法篩選中心基因，表2顯示了分別用這些方法選取的15個(gè)中心基因。最后，把5種方法篩選出的基因進(jìn)行重疊得到了6個(gè)關(guān)鍵基因(圖5)。這6個(gè)基因都存在于這5個(gè)方法中，分別是：Fos、SYK、HCK、SELL、CCR1、NLRP3。

表2 用cytoHubba篩選的15個(gè)中心基因Table 2 The top 15 hub genes rank in cytoHubba

圖5 使用重疊方法對(duì)cytoHubba中篩選出的15個(gè)基因進(jìn)行重疊確定了6個(gè)關(guān)鍵基因Fig.5 Six hub genes were identified by overlapping the first 15 genes in the cytoHubba

2.6 免疫細(xì)胞浸潤分析

利用CIBERSORT算法，研究了PMOP的22個(gè)免疫細(xì)胞集與正常組織免疫滲透的差異。圖6a總結(jié)了3例正常對(duì)照組和10例PMOP患者的結(jié)果。與正常組織相比，PMOP組織通常含有較低比例的活性CD4T記憶細(xì)胞(activated T cells CD4 memory)，見圖6b。

2.7 對(duì)差異基因的驗(yàn)證

NLRP3是我們得到的關(guān)鍵基因，基于GSE7429數(shù)據(jù)集對(duì)NLRP3的表達(dá)進(jìn)行驗(yàn)證(圖7)。這一結(jié)果與我們綜合分析的結(jié)果基本一致。

3 討論

GO分析表明，DEG主要參與生物過程，基因主要在先天免疫反應(yīng)中富集，可能是因?yàn)槿狈S生素D，影響骨骼的代謝所產(chǎn)生的結(jié)果，維生素D受體由大多數(shù)免疫細(xì)胞表達(dá)，包括B和T淋巴細(xì)胞和樹突狀細(xì)胞，其次，免疫細(xì)胞具有活躍的維生素D代謝功能，可以將25(OH)D3局部轉(zhuǎn)移到1,25(OH)2D3中[12]，表明維生素D在調(diào)節(jié)免疫功能方面起著關(guān)鍵作用。此外，之前研究表明先天性免疫系統(tǒng)和腸道微生物群之間的通信障礙可能導(dǎo)致復(fù)雜的疾病[13]，而且腸道微生物群的變化與老年人骨礦物質(zhì)密度的降低有關(guān)[14]，這些研究也側(cè)面的證實(shí)了我們實(shí)驗(yàn)的富集功能結(jié)果，我們推測，這些基因是通過先天免疫系統(tǒng)影響腸道微生物群，進(jìn)而影響雌激素對(duì)破骨細(xì)胞的影響，最終導(dǎo)致骨質(zhì)疏松癥的發(fā)生[15]。

注：a GSE56116數(shù)據(jù)集中13個(gè)樣本中22個(gè)免疫細(xì)胞亞群的相對(duì)百分比；b PMOP與正常對(duì)照組免疫浸潤的差異。正常對(duì)照組為藍(lán)色，PMOP組為紅色，P<0.05。圖6 PMOP患者與正常對(duì)照組的免疫浸潤情況Fig.6 The landscape of immune infiltration between PMOP and normal controls

KEGG通路分析表明，差異基因主要參與甲狀腺素合成等過程。骨骼是T3作用的關(guān)鍵目標(biāo)組織[16]，甲狀腺素短缺和過量都與骨折風(fēng)險(xiǎn)增高有關(guān)[17]。因此，了解T3的分子作用機(jī)制可以指導(dǎo)骨質(zhì)疏松癥的控制和治療。在關(guān)鍵模塊中，我們發(fā)現(xiàn)，該模塊的差異基因在GO富集分析中炎癥反應(yīng)、先天免疫反應(yīng)處富集最為顯著，這與我們之前的功能富集分析結(jié)果一致。此外，我們確定了6個(gè)關(guān)鍵基因：FOS、SYK、HCK、SELL、CCR1、NLRP3。其中NLRP3是最重要的一個(gè)關(guān)鍵基因，此基因在骨侵蝕水平中起著重要作用，NLRP3的功能表達(dá)可能是導(dǎo)致感染部位骨質(zhì)丟失的關(guān)鍵因素[18]。先前的研究發(fā)現(xiàn)將SYK 募集至磷酸化的 ITAM對(duì)破骨細(xì)胞的形成至關(guān)重要[19]，SYK是嗜中性粒細(xì)胞中整合信號(hào)傳導(dǎo)的重要組成部分[20]。HCK和SRC 在成骨細(xì)胞中具有部分重疊的功能，而 HCK 在成骨的表達(dá)中可改善其功能缺陷[21]。CCR1可以改變成骨細(xì)胞和破骨細(xì)胞的功能和分化，以及改變成骨細(xì)胞和脂肪形成之間的不平衡[22]，關(guān)節(jié)炎患者的破骨細(xì)胞前體細(xì)胞高表達(dá) CCR1[23]。因此，對(duì)這些基因的深入了解有助于闡明PMOP形成的關(guān)鍵機(jī)制。然而，我們研究的局限性是，沒有進(jìn)一步的臨床試驗(yàn)來驗(yàn)證它，但它仍然會(huì)對(duì)以后的研究產(chǎn)生很大的幫助。

注：X軸表示基因，Y軸表示基因表達(dá)水平。圖7 驗(yàn)證在GSE7429數(shù)據(jù)P集的關(guān)鍵基因的表達(dá)水平Fig.7 Validation of the expression levels of selected based on GSE7429

綜上，我們確定了308個(gè)差異基因、6個(gè)關(guān)鍵基因和一個(gè)重要網(wǎng)絡(luò)模塊，富集分析了差異基因的信號(hào)通路，通過將可靠的反卷積算法與大規(guī)?；蚪M數(shù)據(jù)相結(jié)合，發(fā)現(xiàn)了PMOP和正常對(duì)照之間的免疫滲透差異。這為臨床診斷和治療提供了可靠的參考價(jià)值。