絕經后骨質疏松生物標志物的篩選及免疫細胞浸潤分析

肖賀方 連笑宇 車凌儀 強碩 鄭稼*

1.鄭州大學人民醫院 河南省人民醫院骨科，河南 鄭州 450003 2.鄭州大學第一附屬醫院婦產科，河南 鄭州 450003

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

1 材料和方法

1.1 數據來源

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

1.2 差異分析

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

1.3 富集分析

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

1.4 構建蛋白質互作網絡(PPI)和篩選關鍵基因

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

1.5 用CIBERSORT 算法分析免疫細胞的浸潤

利用CIBERSORT算法對獲得的數據進行了分析，并得到了22種免疫細胞的比例?；赑值 <0.05 篩選，并計算樣本中每個免疫細胞的百分比。通過R軟件中的vioplot包，分析兩組免疫細胞的不同免疫滲透水平。

1.6 在數據庫中進行驗證

GSE7429 數據集是從GEO數據庫獲得，使用此數據集驗證得到的關鍵基因。

2 結果

2.1 確定差異基因

從數據中篩選出了308個差異基因，其中110個下調，198個上調。圖1和圖2分別顯示了所有DEG的火山圖和前50個下調基因和50個上調基因的熱圖。

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

圖2 紅色表示熱圖中前50上調的差異基因，藍色表示前50個下調的差異基因Fig.2 The heat map of the top 50 upgrade genes in red and 50 downgrade genes in blue

2.2 GO和KEGG的富集分析

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

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

2.3 KEGG通路分析

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

2.4 構建PPI網絡

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

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

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

表1 對關鍵模塊的差異基因進行GO和KEGG富集分析的前10項富集結果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 確定關鍵基因

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

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

圖5 使用重疊方法對cytoHubba中篩選出的15個基因進行重疊確定了6個關鍵基因Fig.5 Six hub genes were identified by overlapping the first 15 genes in the cytoHubba

2.6 免疫細胞浸潤分析

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

2.7 對差異基因的驗證

NLRP3是我們得到的關鍵基因，基于GSE7429數據集對NLRP3的表達進行驗證(圖7)。這一結果與我們綜合分析的結果基本一致。

3 討論

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

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

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

注：X軸表示基因，Y軸表示基因表達水平。圖7 驗證在GSE7429數據P集的關鍵基因的表達水平Fig.7 Validation of the expression levels of selected based on GSE7429

綜上，我們確定了308個差異基因、6個關鍵基因和一個重要網絡模塊，富集分析了差異基因的信號通路，通過將可靠的反卷積算法與大規模基因組數據相結合，發現了PMOP和正常對照之間的免疫滲透差異。這為臨床診斷和治療提供了可靠的參考價值。