女性骨關(guān)節(jié)炎患者外周血單個(gè)核細(xì)胞基因表達(dá)譜的生物信息學(xué)分析

楊 威，袁普衛(wèi)，杜龍龍，韓清民

（1. 廣州中醫(yī)藥大學(xué)第三臨床醫(yī)學(xué)院，廣東廣州 510405；2. 陜西中醫(yī)藥大學(xué)，陜西咸陽(yáng) 712046；3. 廣州中醫(yī)藥大學(xué)第三附屬醫(yī)院，廣東廣州 510405）

骨關(guān)節(jié)炎(osteoarthritis, OA)是一種以軟骨不可逆性退變?yōu)橹鞯淖畛Ｒ?jiàn)關(guān)節(jié)病，使受累關(guān)節(jié)負(fù)荷疼痛和功能受限，好發(fā)于中老年女性，嚴(yán)重影響患者生活。OA 主要病理特征為關(guān)節(jié)軟骨降解、滑膜炎癥及其周圍組織重塑[1]，與遺傳、衰老及慢性炎癥等多因素有關(guān)，但其病因病機(jī)不明。隨著人口老齡化及肥胖加劇，患病率持續(xù)增長(zhǎng)已成全球共識(shí)。目前，臨床沒(méi)有監(jiān)測(cè)OA 發(fā)生或進(jìn)展的特異靶點(diǎn)，不僅缺乏對(duì)OA 的早期診斷，也無(wú)減緩疾病進(jìn)展的方法。因此，探索對(duì)OA 敏感且實(shí)用的分子標(biāo)志物十分迫切。而血液作為臨床信息的重要來(lái)源，提供包括代謝物、蛋白質(zhì)和核酸等信息，隨著分子生物學(xué)的發(fā)展，與OA 相關(guān)的血液生物學(xué)指標(biāo)也逐漸增多，尤其是關(guān)注外周血單個(gè)核細(xì)胞(peripheral blood mononuclear cells, PBMCs)的生物標(biāo)志物。已有國(guó)內(nèi)外研究鑒別出OA 患者PBMCs基因表達(dá)譜中可能用于診治的差異基因[1-4]，但都是對(duì)所有OA 患者PBMCs 基因分析，沒(méi)有區(qū)分性別，也沒(méi)有采用基因集合富集分析(Gene Set enrichment analysis, GSEA)，而且所有分析結(jié)果都不一致。

本研究分析數(shù)據(jù)集中所有女性外周血樣本，基于基因集富集分析(gene set enrichment analysis,GSEA)尋找與OA相關(guān)的目的基因，去重后進(jìn)行經(jīng)典的基因本論(gene ontology,GO)富集、京都基因與基因組百科全書(shū)(Kyoto Encyclopedia of Genes and Genomes, KEGG)通路富集和蛋白互作可視化網(wǎng)絡(luò)(protein-protein interaction network,PPI network)分析，運(yùn)用MCODE、cytohubba 和centiscape 插件篩選核心基因，為尋找OA 診治的適宜分子靶點(diǎn)提供理論基礎(chǔ)。

1 資料與方法

1.1 數(shù)據(jù)來(lái)源

從GEO（https://www.ncbi.nlm.nih.gov/geo)數(shù)據(jù)庫(kù)查找OA 患者PBMCs 基因表達(dá)芯片GSE48556[2],包括106 例OA 患者血液樣本和33 例正常志愿者（NC）樣本，由于僅1 例男性O(shè)A 患者，樣本量過(guò)少，無(wú)法進(jìn)行性別對(duì)比，故予剔除。本研究針對(duì)數(shù)據(jù)集中所有女性患者，納入情況見(jiàn)表1，所處平臺(tái)為GPL6947 Illumina HumanHT-12 V3.0 expression beadchip，類型為表達(dá)譜芯片。本研究共有48 802 條數(shù)據(jù)，25 159個(gè)基因。

表1 納入樣本情況Tab.1 Characteristics of the included samples

1.2 數(shù)據(jù)質(zhì)量與差異分析

用GEO2R（http://ncbi.nlm.nih.gov/geo/geo2r）對(duì)數(shù)據(jù)集進(jìn)行分組及分析，使用R 語(yǔ)言的limma 包對(duì)樣本進(jìn)行一致性檢驗(yàn)和主成分分析(principal component analysis, PCA)，用以檢測(cè)樣本的均一性、可重復(fù)性和組間差異性。

1.3 GSEA 分析

將原始數(shù)據(jù)經(jīng)R 語(yǔ)言計(jì)算出相應(yīng)的P.adjust、qvalue、Pvalue 和log2基因表達(dá)差異倍數(shù)(fold-change,FC)等，用clusterProfiler 包進(jìn)行GSEA 分析，于分子特征數(shù)據(jù)庫(kù)官網(wǎng)(http://www.gsea-msigdb.org/gsea/msigdb/index.jsp)分別選擇基因集C2.CP.KEGG.v7.4、C5.GO.BP/CC/MF.v7.4做相應(yīng)分析。用GSEA分析找到富集的生物學(xué)過(guò)程或通路，將對(duì)應(yīng)的所有l(wèi)eading-edge 基因累加后去重。

1.4 目的基因篩選

將得到的leading-edge 基因子集用excel 2010 在經(jīng)R 語(yǔ)言預(yù)處理過(guò)的原始數(shù)據(jù)中找到相對(duì)應(yīng)的P.adjust、qvalue、Pvalue 和log2FC 值，篩 選 條 件 是P.adjust＜0.05&log2FC|＞0.2[5]，共得到目的基因292 個(gè)。

1.5 GO 和KEGG 富 集 分 析

使用DAVID 數(shù)據(jù)庫(kù)(https://david.ncifcrf.gov/tools.jsp）對(duì)篩選出的292 個(gè)基因進(jìn)行GO 和KEGG通路富集。GO 分析包括生物學(xué)過(guò)程(biological process, BP)、細(xì)胞成分(cellular component, CC)和分子功能(molecular function, MF)，用于預(yù)測(cè)蛋白質(zhì)的功能。KEGG 通路分析出分子間相互作用、反應(yīng)和構(gòu)建關(guān)系網(wǎng)絡(luò)。P＜0.05 認(rèn)為差異具有統(tǒng)計(jì)學(xué)意義。

1.6 PPI 網(wǎng)絡(luò)構(gòu)建及模塊分析

將目的基因輸入STRING(https://string-db.org/)在線平臺(tái)，設(shè)置combined score＞0.4，然后將計(jì)算結(jié)果導(dǎo)入Cytoscape-v3.8.2 中，構(gòu)建PPI 網(wǎng)絡(luò)。通過(guò)MCODE 插件分析出PPI 網(wǎng)絡(luò)中連接最為緊密的基因集簇，使用cytohubba 插件采用MCC[3]和Degree算法篩選出連結(jié)度最高的前10 個(gè)基因，再聯(lián)合centiscape 插件中Degree、Betweenness、Closeness 三個(gè)參數(shù)以及客觀的P值、P.adjust 和log2FC 值綜合分析篩選核心基因。

2 結(jié) 果

2.1 數(shù)據(jù)質(zhì)量與差異分析

將GSE48556 數(shù)據(jù)集所有女性分為OA 組和NC組，分別有105、24 例，標(biāo)準(zhǔn)化后各樣本中位數(shù)基本在一條水平線上（圖1）。二維及三維PCA 分析，分組間差異不夠明顯，說(shuō)明PCA 對(duì)組間差異解釋有限（圖2）。ggplot2 包制作全基因組火山圖，F(xiàn)C 值設(shè)定為1.2，P＜0.05 有統(tǒng)計(jì)學(xué)意義，顯著上調(diào)有326 個(gè)基因，顯著下調(diào)有385 個(gè)基因，ComplexHeatmap 包制作熱圖，顯示表達(dá)譜中高、低表達(dá)各top20 基因（圖3）。

圖3 GSE48556 數(shù)據(jù)集差異表達(dá)的火山圖和熱圖Fig.3 The differentially expressed volcano map and heat map of the GSE48556 data set

A、B 分別為歸一化前后，綠色為NC 組代表正常對(duì)照，藍(lán)色為OA 組。

A、B 圖分別為二維或三維，PC1、PC2、PC3 分別表示主成分1、2、3，NC 為正常對(duì)照組。

2.2 GSEA 分析

GSEA 分析找到富集的生物學(xué)過(guò)程或通路，篩選條件為：BP 和CC 以P.adjust＜0.05 &qvalue＜0.25，MF 和KEGG 以P＜0.05 &qvalue＜0.25，分別篩到80、37、22、16 條，將對(duì)應(yīng)的所有l(wèi)eading-edge 基因累加后去重，一共得到3 120 個(gè)基因。富集到的BP 多與固有免疫和適應(yīng)性免疫調(diào)節(jié)等有關(guān)；CC 多為細(xì)胞核仁、核膜、囊泡和內(nèi)質(zhì)網(wǎng)等；MF 多與酶、轉(zhuǎn)錄因子或組蛋白等的結(jié)合有關(guān)；KEGG 主要與移植物抗宿主病、抗原呈遞和自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性等有關(guān)，以下分別展示與OA 相關(guān)的功能或通路（圖4）。

圖4 基因芯片GSE48556 數(shù)據(jù)集GSEA 富集分析Fig.4 GSEA enrichment analysis of gene chip GSE48556 data set

2.3 目的基因

經(jīng)過(guò)篩選得到的leading-edge 基因子集，有3 120個(gè)基因，其中上調(diào)基因1 195 個(gè)，下調(diào)基因1 925 個(gè)；考慮到功能基因的差異和顯著性，滿足P.adjust＜0.05 & |log2FC|＞0.2 的有292 個(gè)目的基因，其中上調(diào)81 個(gè)，下調(diào)211 個(gè)，分別做熱圖對(duì)比展示（圖5）。

圖5 leading-edge 基因子集和目的基因熱圖Fig.5 The heat maps of leading edge gene subset and target genes

2.4 GO 功能和KEGG 通路富集分析

對(duì)目的基因進(jìn)行GO 和KEGG 富集分析，共富集到289 條，其中在滿足P.adjust＜0.05 &qvalue＜0.2條件下，BP 共有124 條，CC 共有12 條，MF 共有8 條，KEGG 共有11 條。BP 主要富集在NIK/NF-κB 信號(hào)的調(diào)節(jié)、單核細(xì)胞增殖、凋亡信號(hào)通路的調(diào)控、TNF介導(dǎo)的信號(hào)通路和Wnt 信號(hào)通路的調(diào)控等；CC 富集在細(xì)胞核、核膜和胞質(zhì)連接等；MF 富集在蛋白質(zhì)結(jié)合、分子接頭、細(xì)胞趨化因子受體活性或調(diào)節(jié)等；KEGG 富集在自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性、TNF信號(hào)通路、MAPK 信號(hào)通路、細(xì)胞凋亡、T 細(xì)胞受體信號(hào)通路、FoxO 信號(hào)通路、NF-κB 信號(hào)通路和趨化因子信號(hào)通路等。將P.adjust 的值從小到大排序，結(jié)合與OA 的相關(guān)性展示BP 和KEGG，分別有12 條和10 條，CC 和MF 全部展示，見(jiàn)表2、圖6。

圖6 GO 和KEGG 富集分析氣泡圖Fig.6 Bubble charts of GO and KEGG enrichment analyses

表2 GO 和KEGG 富集分析Tab.2 GO and KEGG enrichment analyses

2.5 PPI 網(wǎng)絡(luò)分析及模塊分析結(jié)果

PPI 網(wǎng)絡(luò)圖共有235 個(gè)節(jié)點(diǎn)，687 條互作關(guān)系，經(jīng)MCODE 篩選出7 個(gè)集簇，展示得分較高的前4 個(gè)集簇，其中連接度最高的基因集簇，包括11 個(gè)節(jié)點(diǎn)，47條互作關(guān)系（圖7）。然后使用cytohubba 插件采用MCC 和Degree 算法篩選出連結(jié)度最高的前10 個(gè)基因，再聯(lián)合centiscape 插件中度(degree)、中介中心性(betweenness)和緊密度(closeness)等參數(shù)以及對(duì)應(yīng)的P.adjust、Pvalue 和log2FC 值篩選核心基因，最后結(jié)合與OA 緊密相關(guān)的分子生物學(xué)意義綜合分析確定8 個(gè)核心基因?yàn)椋篗APK1、IL10、PTGS2、IL18、GSK3B、NFKBIA、TNFRSF1A、EGR1（圖8 和表3）。

表3 核心基因關(guān)系網(wǎng)絡(luò)的拓?fù)鋵W(xué)分析Tab.3 Topological analysis of core gene relationship network

圖7 PPI 網(wǎng)絡(luò)及模塊分析圖Fig.7 PPI network and module analysis diagram

3 討 論

OA 是中老年女性患者疼痛和致殘的主要原因之一，不僅影響正常生活，也嚴(yán)重限制社交活動(dòng)。目前以姑息治療為主，尚無(wú)阻斷或逆轉(zhuǎn)關(guān)節(jié)退變的方法。隨著疾病的不可逆性進(jìn)展，多數(shù)患者需要接受終末期手術(shù)治療。因此，及時(shí)診治尤為重要。而血液作為傳統(tǒng)的臨床檢測(cè)樣本，因其近乎無(wú)創(chuàng)和易獲得性，得到研究者的青睞。血液中PBMCs，是外周循環(huán)血中具有單個(gè)核細(xì)胞的總稱，單個(gè)核是相對(duì)多個(gè)核或分葉核粒細(xì)胞而言，除間充質(zhì)干細(xì)胞、造血干細(xì)胞和多潛能細(xì)胞等少數(shù)祖細(xì)胞群體外[6]，主要包括淋巴細(xì)胞、單核細(xì)胞(Monocyte,Mo)、樹(shù)突狀細(xì)胞和NK 細(xì)胞等，Mo 可以發(fā)育為各種組織巨噬細(xì)胞，PBMCs為免疫系統(tǒng)的主要組成部分，參與人體免疫應(yīng)答。近年來(lái)，人們開(kāi)始關(guān)注OA 外周血成分的變化。PONCHEL等[7]發(fā)現(xiàn)，與健康人相比，OA 患者血中CD4+T 細(xì)胞和B 細(xì)胞更少，而CD8+T 細(xì)胞和記憶細(xì)胞更多，并且所表現(xiàn)的免疫功能障礙與衰老無(wú)關(guān)。臨床上，王猛等[8]發(fā)現(xiàn)膝OA 患者外周血中性粒細(xì)胞-淋巴細(xì)胞比率和單核細(xì)胞-淋巴細(xì)胞比率明顯高于健康志愿者，且對(duì)全膝置換術(shù)后關(guān)節(jié)功能評(píng)估有參考意義。

本研究基于GSE48556 數(shù)據(jù)集中所有女性血液樣本，105 例髖或膝OA 患者和24 例正常對(duì)照，所有樣本年齡和體質(zhì)量指數(shù)基線指標(biāo)無(wú)差異，二維和三維PCA 分析顯示樣本間距離小，可能與遺傳背景接近有關(guān)，絕大部分OA 患者有疼痛，約一半服用非甾體類抗炎藥。傳統(tǒng)Log2FC 閾值篩選差異基因法沒(méi)有考慮差異微弱的基因，可能會(huì)造成很大的偏倚，尤其對(duì)血中差異基因表達(dá)微弱疾病的鑒別。而GSEA 富集法是對(duì)所有表達(dá)水平的基因分析，將其與預(yù)定義的基因集比較，計(jì)算每個(gè)基因的富集分?jǐn)?shù)，不僅可確定待測(cè)基因組在特定基因集中的表達(dá)狀況及是否存在統(tǒng)計(jì)學(xué)差異，還可篩選出對(duì)富集得分貢獻(xiàn)最大的基因子集。經(jīng)過(guò)GSEA 富集篩選得到的3 120 個(gè)功能基因二次分析，發(fā)現(xiàn)很多功能基因沒(méi)有統(tǒng)計(jì)學(xué)意義，log2FC值太小，以P.adjust＜0.05 &|log2FC|＞0.2 再次篩選，共得到目的基因292 個(gè)行GO、KEGG 和PPI 分析。

GO 分析顯示，目的基因主要富集在“泛素依賴性蛋白質(zhì)分解代謝過(guò)程的調(diào)節(jié)”、“NIK/NF-κB 的調(diào)節(jié)”、“單核細(xì)胞增殖”、“凋亡信號(hào)通路的調(diào)控”、“細(xì)胞對(duì)機(jī)械刺激的反應(yīng)”、“TNF 介導(dǎo)的信號(hào)通路”、“Wnt信號(hào)通路的調(diào)控”、“MAP 激酶活性的調(diào)節(jié)”和“自噬的調(diào)節(jié)”等生物功能上。已有研究報(bào)道過(guò)OA 類似的病理機(jī)制。如LIU 等[9]認(rèn)為泛素和泛素樣修飾酶促機(jī)制在軟骨細(xì)胞穩(wěn)態(tài)方面發(fā)揮重要作用，包括調(diào)節(jié)細(xì)胞信號(hào)轉(zhuǎn)導(dǎo)和蛋白質(zhì)穩(wěn)定性，處理細(xì)胞應(yīng)激和炎癥，以及維持軟骨細(xì)胞分化和存活，可能成為治療OA 新的靶點(diǎn)。王猛等[8]發(fā)現(xiàn)OA 患者外周血Mo 明顯高于健康志愿者，說(shuō)明單核細(xì)胞增殖與OA 有關(guān)，并且，外周血Mo 是滑膜巨噬細(xì)胞的來(lái)源細(xì)胞，而后者可能是OA 炎 癥 的 發(fā) 起 者[10]。TER HEEGDE 等[11]研 究 發(fā) 現(xiàn)，機(jī)械負(fù)荷導(dǎo)致細(xì)胞應(yīng)激反應(yīng)和局部炎癥，所誘發(fā)的疼痛與軟骨損害的OA 樣改變有關(guān)。NF-κB、凋亡信號(hào)通路、TNF 相關(guān)通路、Wnt 信號(hào)通路、MAP 激酶和自噬調(diào)節(jié)是常見(jiàn)的OA 病理機(jī)制，多與軟骨細(xì)胞凋亡和滑膜細(xì)胞炎癥有關(guān)，其中，NF-κB、TNF 相關(guān)通路和自噬與炎癥機(jī)制密切相關(guān)，凋亡信號(hào)通路、Wnt 信號(hào)通路和MAP 激酶調(diào)節(jié)與軟骨細(xì)胞凋亡更為緊密，兩者之間又相互影響。KEGG 主要富集于以下4 條與OA 相關(guān)的通路：自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性、TNF 信號(hào)通路、MAPK 信號(hào)通路和細(xì)胞凋亡。而關(guān)于TNF 信號(hào)通路、MAPK 信號(hào)通路和細(xì)胞凋亡在OA 病機(jī)中的研究較多，且在很多文獻(xiàn)都有報(bào)道。關(guān)于“自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性”，F(xiàn)ENG 等[1]對(duì)GSE48556 數(shù)據(jù)集分析時(shí)也發(fā)現(xiàn)其是OA 差異基因主要富集通路之一，但目前尚未發(fā)現(xiàn)其參與OA 病機(jī)的相關(guān)實(shí)驗(yàn)研究。KEGG 富集分析進(jìn)一步說(shuō)明炎癥和凋亡與OA 的病理變化密切相關(guān)。

通過(guò)PPI網(wǎng)絡(luò)與centiscape插件篩選出8個(gè)核心基因：MAPK1、IL10、PTGS2、IL18、GSK3B、NFKBIA、TNFRSF1A、EGR1，其中僅MAPK1 和TNFRSF1A上調(diào)，其余基因都是下調(diào)。EGR1、PTGS2和IL18 都被之前的分析篩選為核心基因[1,4]。MAPK1 是MAP激酶家族一員，MAPK1 磷酸化加重軟骨細(xì)胞肥大分化，在炎癥因子共同作用下引起軟骨細(xì)胞凋亡。IL10主要由單核細(xì)胞產(chǎn)生，具有強(qiáng)大的抗炎功能，限制炎癥引起的過(guò)度組織破壞。CAMERON 等[12]在評(píng)估骨髓源性間充質(zhì)干細(xì)胞(BM-MSC)和IL-10 在OA 中的組合作用的細(xì)胞實(shí)驗(yàn)中發(fā)現(xiàn)，IL-10 過(guò)表達(dá)可能會(huì)增強(qiáng)BM-MSC 介導(dǎo)的T 細(xì)胞抑制，但沒(méi)有發(fā)現(xiàn)腺相關(guān)病毒-IL10 轉(zhuǎn)導(dǎo)的BM-MSC 中炎性驅(qū)動(dòng)軟骨降解的顯著調(diào)節(jié)。GSK3B 是糖原合酶激酶-3β，葡萄糖穩(wěn)態(tài)的負(fù)調(diào)節(jié)劑，參與能量代謝、炎癥、內(nèi)質(zhì)網(wǎng)應(yīng)激和細(xì)胞凋亡等，HUANG 等[13]之前的研究中已將其篩選為核心基因，通過(guò)抑制WNT/β-catenin 通路維持軟骨表型和軟骨細(xì)胞外基質(zhì)，拮抗軟骨細(xì)胞凋亡。NFKBIA是NF-κB 抑制劑α，抑制參與炎癥反應(yīng)的NF-κB 復(fù)合物，TNFRSF1A 是腫瘤壞死因子受體超家族成員1A，與TNFα 受 體 結(jié) 合 抑 制 炎 癥，XU 等[14]發(fā) 現(xiàn)，TNF-α 刺激下的軟骨細(xì)胞NF-κB 通路高度富集，木蘭素可以抑制NFKBIA 和TNFRSF1A 的表達(dá)而抑制OA 中軟骨細(xì)胞的炎性反應(yīng)。

然而，結(jié)合臨床意義分析核心基因時(shí)發(fā)現(xiàn)：MAPK1、EGR1 和GSK3B 主要與軟骨細(xì)胞凋亡有關(guān)，前兩者正相關(guān)，后者負(fù)相關(guān)，但是正相關(guān)的EGR1表達(dá)量下調(diào)，與OA 軟骨細(xì)胞凋亡增加不符，可能與機(jī)體自發(fā)的負(fù)反饋調(diào)節(jié)有關(guān)；IL18、PTGS2、IL10、NFKBIA 和TNFRSF1A 主要與炎癥反應(yīng)有關(guān)，前兩者為促炎因子，后三者為抗炎因子，但是IL18 和PTGS2 的表達(dá)量都下調(diào)，可能與患者使用非甾體類抗炎藥有關(guān)，也可能與機(jī)體自發(fā)保護(hù)反應(yīng)有關(guān)。另外，本研究有以下不足：①對(duì)數(shù)據(jù)集的檢索可能不完全，沒(méi)有納入全部相關(guān)資料；②沒(méi)有建立合適的數(shù)學(xué)模型，仍然存在人為設(shè)置閾值，手工挑選目的基因；③沒(méi)有進(jìn)一步的實(shí)驗(yàn)驗(yàn)證。總之，本研究發(fā)現(xiàn)：①自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性、TNF 信號(hào)通路、MAPK 信號(hào)通路、細(xì)胞凋亡等生物學(xué)過(guò)程與OA 病機(jī)密切相關(guān)；②MAPK1、IL10、PTGS2、IL18、GSK3B、NFKBIA、TNFRSF1A、EGR1 八個(gè)核心基因與OA炎癥和細(xì)胞凋亡高度相關(guān)，可能是PBMCs 中區(qū)別于健康人的有效生物標(biāo)志物，有望通過(guò)外周血檢測(cè)識(shí)別OA。