轉(zhuǎn)錄組分析鑒定腰椎椎間盤退行性變中潛在免疫相關(guān)生物標志物

李志超，王 磊，薛景才，王文波，陳 煜，齊軍強，許 鵬

1.山東中醫(yī)藥大學第二附屬醫(yī)院骨科，濟南 250002

2.山東中醫(yī)藥大學研究生院，濟南 250001

3.海軍軍醫(yī)大學研究生院，上海 200433

4.海軍軍醫(yī)大學長征醫(yī)院骨科，上海 200003

以椎間盤退行性變（IDD）為特征的椎間盤疾病，給患者帶來了生理和心理上的痛苦，給家庭和社會造成巨大的經(jīng)濟負擔［1］。IDD 的致病因素包括營養(yǎng)物質(zhì)減少、基質(zhì)降解、細胞凋亡、炎性因子增加、免疫應答、遺傳因素和機械負荷等［2］，但其具體機制目前仍不清楚。IDD 與免疫之間有很強的相關(guān)性［3］，椎間盤的駐留細胞（纖維環(huán)和髓核）和浸潤的免疫細胞分泌的炎性細胞因子是引發(fā)IDD 的關(guān)鍵因素［4］。發(fā)生退行性變的椎間盤使巨噬細胞傾向于促炎性極化，這似乎加劇了IDD［5］。

競爭性內(nèi)源RNA（ceRNA）是一種新的基因表達調(diào)控模式。在調(diào)控過程中，一些長鏈非編碼RNA（lncRNA）和微RNA（miRNA）結(jié)合，lncRNA 起到miRNA 海綿的作用，降低miRNA 對其靶基因的抑制作用，使靶基因在細胞中的表達量增加［6］。有研究［7-8］表明，IDD 的病理變化與lncRNA、miRNA 及其靶點的失衡以及許多不同細胞生物過程有關(guān)。失衡的ceRNA 相互作用軸可能是治療IDD 的關(guān)鍵目標［9］。本研究通過加權(quán)基因共表達分析（WGCNA）構(gòu) 建IDD 中l(wèi)ncRNA-miRNA-mRNA 的ceRNA 網(wǎng) 絡，篩選與IDD 免疫相關(guān)的關(guān)鍵基因，為尋找IDD 的新靶點和治療策略提供參考。

1 材料和方法

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

從基因表達匯編（GEO）數(shù)據(jù)庫（https：//www.ncbi.nlm.nih.gov/geo）GSE67567 數(shù) 據(jù) 集 下 載IDD 人類退行性髓核組織和對照的髓核組織微陣列芯片數(shù)據(jù)，GSE67567 數(shù)據(jù)集包含GSE56081（lncRNA 和mRNA 數(shù)據(jù)）和GSE63492（miRNA）2 個子集，各包括來自5名健康志愿者和5例IDD患者的髓核樣本。

1.2 數(shù)據(jù)預處理和差異表達分析

使用GEO 的在線分析工具GEO2R 獲得標準化的差異表達基因。根據(jù)GPL15314、GPL19449文件和GEO中的注釋信息，將探針轉(zhuǎn)換為相應的基因符號。對于差異表達的miRNA，將閾值設置為|log2FC| ＞ 2、調(diào)整P＜ 0.05；對于差異表達的lncRNA，將閾值設置為|log2FC| ＞ 1、調(diào)整P＜ 0.05。

1.3 共表達模塊建立和基因集富集分析（GSEA）

使用Image GP工具（https：//www.ehbio.com/Cloud Platform/front/#/）分 析GEO 中的mRNA 表達矩陣。首先構(gòu)建差異表達基因的共表達網(wǎng)絡，將其轉(zhuǎn)換為拓撲重疊矩陣，通過分層聚類和動態(tài)樹切割算法，使用模塊特征基因和模塊關(guān)系來確定與表型相關(guān)的共表達模塊，最小模塊大小為30［10］。使用Enrichr數(shù)據(jù)庫（https：//maayanlab.cloud/Enrichr/）分析每個模塊中的mRNA，進行京都基因與基因組百科全書（KEGG）通路和基因本體（GO）生物過程的GSEA，篩選出與免疫調(diào)節(jié)關(guān)系最為密切的模塊。

1.4 ceRNA 網(wǎng)絡建立和GSEA

通過Targetscan 數(shù)據(jù)庫（https：//www.targetscan.org）和miRTarbase 數(shù)據(jù)庫（http：//mirtarbase.mbc.nctu.edu.tw/php/index.php）篩 選miRNA 靶 向 的mRNA，將這些mRNA 與上述所選模塊中的差異表達mRNA 取交集獲得免疫相關(guān)mRNA。同理，通過StarBase 數(shù)據(jù)庫（http：//starbase.sysu.edu.cn/index.php）預測差異表達miRNA 調(diào)控的lncRNA。最后，使用Cytoscape 3.7.2建立lncRNA-miRNA-mRNA ceRNA網(wǎng) 絡。使 用Enrichr 對ceRNA 網(wǎng) 絡 中 的mRNA 的KEGG 通路和GO 生物過程進行GSEA。

1.5 免疫浸潤分析

使用ImmuCellAI（http：//bioinfo.life.hust.edu.cn/web/ImmuCellAI/）［11］比較GEO 數(shù)據(jù)庫中IDD 組和對照組的免疫浸潤特征。ImmuCellAI 是一種基于轉(zhuǎn)錄組測序或芯片工具估計免疫細胞浸潤豐度數(shù)據(jù)的工具，可以判斷18 種T 細胞和6 種其他類型免疫細胞的比例。

1.6 與CD4 幼稚T 細胞（CD4 naive）相關(guān)的ceRNA子網(wǎng)絡的構(gòu)建

使用Pearson 相關(guān)分析確定ceRNA 網(wǎng)絡中CD4 naive 浸潤與mRNA 和lncRNA 表達的相關(guān)性，以相關(guān)系數(shù)絕對值＞ 0.8、P＜ 0.01 為條件篩選mRNA和lncRNA［12］。同法分析mRNA 和lncRNA 之間的相關(guān)性，篩選出關(guān)鍵lncRNA。利用篩選出的mRNA 和lncRNA，在之前的ceRNA網(wǎng)絡基礎(chǔ)上構(gòu)建與CD4 naive相關(guān)的ceRNA 子網(wǎng)絡。

1.7 確定關(guān)鍵靶點的通路

使用關(guān)鍵lncRNA 表達量的中位數(shù)對所有表達量進行分組，并使用進行GO 生物過程和KEGG 通路的GSEA。通過獲得的標準化富集評分（NES）及其相應的P值（NOMP-val）篩選關(guān)鍵的生物過程和途徑。

1.8 驗證

使用GEO 數(shù)據(jù)庫中的另一個IDD 數(shù)據(jù)集GSE124272（包括8 例IDD 患者和8 名健康志愿者），通過IDD 和對照組基因表達的比較，分析差異表達基因，對獲得的關(guān)鍵基因進行驗證。

2 結(jié) 果

2.1 IDD 中mRNA、lncRNA 和miRNA的差異表達

通過GEO2R 分析GSE56081 和GSE63492 數(shù)據(jù)集中的數(shù)據(jù)，經(jīng)過校正，得到15 220 個mRNA、370個miRNA 和623 個與IDD 相關(guān)的lncRNA。進一步根據(jù)|log2FC| ＞ 2、調(diào)整P＜ 0.05 篩選出35 個差異表達miRNA，均為低表達；根據(jù)|log2FC| ＞ 1、調(diào)整P＜ 0.05篩選出87個差異表達lncRNA，其中45個為低表達，42 個為高表達（圖1）。

圖1 差異表達基因Fig.1 Differentially expressed genes

2.2 共表達模塊構(gòu)建及GSEA

通過WGCNA 完成IDD 和對照組的共表達模塊構(gòu)建，通過聚類分析和動態(tài)切樹算法，得到5 個不同顏色的模塊，并識別出與臨床特征相關(guān)的共表達模塊（圖2）。

圖2 WGCNA 結(jié)果Fig.2 Results of WGCNA

對每個模塊的KEGG 通路和GO 生物過程進行GSEA，結(jié)果表明，與免疫調(diào)節(jié)最相關(guān)的是藍綠色（turquoise）模塊，總共包含1 808 個mRNA。藍綠色模塊的KEGG 通路主要包括剪接體、溶酶體、RNA降解、亨廷頓病、人T 細胞白血病病毒1 型感染、氧化磷酸化等，GO 生物學過程主要包括參與免疫反應的中性粒細胞活化、纖溶酶原激活的調(diào)控、巨自噬調(diào)節(jié)、TORC1 信號的調(diào)節(jié)、中性粒細胞介導的免疫、中性粒細胞脫顆粒等（圖3）。

圖3 藍綠色模塊的KEGG通路和GO生物過程Fig.3 KEGG pathways and GO biological processes of turquoise module

2.3 ceRNA 網(wǎng)絡構(gòu)建和GSEA

使用Targetscan 和miRTarbase 數(shù)據(jù)庫預測受差異表達miRNA 調(diào)節(jié)的mRNA，Targetscan 數(shù)據(jù)庫中獲得2 617 個、miRTarbase 數(shù)據(jù)庫中獲得3 940 個mRNA，取交集后獲得160 個mRNA，再與免疫相關(guān)模塊中的mRNA 取交集得到25 個mRNA。使用StarBase 數(shù)據(jù)庫預測出差異表達miRNA 調(diào)節(jié)的797個lncRNA，將其與差異表達lncRNA 取交集后獲得8 個lncRNA。使 用Cytoscape 3.7.2 構(gòu) 建lncRNAmiRNA-mRNA ceRNA 網(wǎng)絡（圖4a），獲得ceRNA 網(wǎng)絡中的mRNA，并對其KEGG 通路和GO 生物過程進行GSEA。結(jié)果表明，GO 生物過程主要包括T 細胞耐受誘導調(diào)控、細胞對粒細胞-巨噬細胞集落刺激因子刺激的反應、CD4+α-β T細胞活化的正調(diào)節(jié)、細胞對轉(zhuǎn)化的反應生長因子刺激的反應、α-β T細胞分化的正調(diào)節(jié)等（圖4b），KEGG通路主要包括Fox O信號通路、細胞衰老、SNARE 在囊泡轉(zhuǎn)運中的相互作用、MAPK 信號通路（圖4c）。

圖4 免疫相關(guān)ceRNA 網(wǎng)絡及GSEAFig.4 Immune-related ceRNA network and GSEA

2.4 與CD4 幼稚T 細胞相關(guān)的ceRNA 子網(wǎng)絡構(gòu)建及GSEA

使用ImmuCellAI 進行免疫細胞浸潤分析，結(jié)果顯示，CD4 naive、CD8 幼稚T 細胞（CD8 naive）、衰竭T 細胞（Tex）、自然調(diào)節(jié)性T 細胞（nTreg）、誘導調(diào)節(jié)性T 細胞（iTreg）、輔助性T 細胞1（Th1）、輔助性T 細胞17（Th17）、濾泡輔助性T 細胞（Tfh）、中央記憶細胞（Tcm）、自然殺傷細胞（NKT）、黏膜相關(guān)不變T 細胞（MAIT）、樹突狀細胞（DC）、單核細胞（monocyte）、γ-σ T 細胞（Tgd）、CD4 T 細胞在IDD組與對照組間差異有統(tǒng)計學意義（P＜ 0.05，圖5a），其中以CD4 naive T 細胞差異最為顯著，IDD 組豐度顯著低于對照組（P＜ 0.05）。

圖5 免疫浸潤分析及l(fā)ncRNA 與mRNA 的Pearson 相關(guān)分析Fig.5 ImmuCellAI immune infiltration analysis and Pearson correlation analysis of lncRNAs and mRNAs

利 用ceRNA 網(wǎng) 絡 中mRNA 和lncRNA 的 表 達量篩選與CD4 naive T 細胞免疫浸潤相關(guān)的關(guān)鍵基因，經(jīng)Pearson 相關(guān)分析，發(fā)現(xiàn)8 個lncRNA 和25 個mRNA 與CD4 naive T細胞浸潤具有強相關(guān)性，其中1 個lncRNA（LINC00641）與10 個mRNA 呈強 正相關(guān)性（圖5b）。

在ceRNA網(wǎng)絡基礎(chǔ)上構(gòu)建一個與CD4 naive T細胞相關(guān)的ceRNA子網(wǎng)絡（圖6a），關(guān)鍵mRNA可能是UHMK1、ZFP36L2、ZCCHC3、ZBTB20。在GSE56081數(shù)據(jù)集的GSEA 中，使用LINC00641 的中位表達值進行分組，明確關(guān)鍵因子的調(diào)控途徑。結(jié)果顯示，許多機制與免疫密切相關(guān)，這些相關(guān)的KEGG 通路包括嗅覺轉(zhuǎn)導、轉(zhuǎn)化生長因子信號通路、溶酶體、胞質(zhì)DNA 感應通路、剪接體等（圖6b）；GO 生物過程主要包括負調(diào)控細胞基質(zhì)黏附、白細胞遷移、細胞對細胞因子刺激的反應、炎性反應、白細胞遷移的調(diào)節(jié)等（圖6c）。

圖6 與CD4 naive相關(guān)的ceRNA 子網(wǎng)絡及GSEAFig.6 ceRNA sub-network associated with CD4 naive and GSEA

2.6 驗證結(jié)果

通過GEO2R 分析GSE124272 數(shù)據(jù)集的差異表達基因，結(jié)果顯示，LINC00641、UHMK1、ZFP36L2、ZCCHC3和ZBTB20這5個關(guān)鍵基因在IDD組和對照組中均存在差異表達（表1）。

表1 GSE124272數(shù)據(jù)集中關(guān)鍵基因的差異表達Tab.1 Differential expression of key genes in GSE124272 dataset

3 討 論

IDD 是臨床常見病，發(fā)生機制復雜［13］。一項薈萃分析［14］表明，突出的椎間盤會自發(fā)重吸收，即IDD 患者未經(jīng)手術(shù)治療或侵入性治療，脫出的髓核自發(fā)消失或明顯縮小，但現(xiàn)有文獻報道多為病例報告或回顧性研究，其機制可能與炎性反應和免疫調(diào)節(jié)有關(guān)［15］。因此，尋找新的免疫靶點和新的治療方法具有重要意義。近年來，基因表達譜的研究為越來越多的疾病篩選新型生物標志物帶來了新思路與新方法［16］，本研究旨在探索IDD 與免疫細胞的相關(guān)性并篩選相關(guān)的特征基因。

本研究通過基因表達譜鑒定差異表達基因，使用WGCNA 和GSEA 篩選免疫相關(guān)模塊，并構(gòu)建了 由8 個lncRNA、21 個miRNA 和25 個mRNA 組成的ceRNA 網(wǎng)絡。mRNA 的GO 生物過程與KEGG通路分析表明其與免疫關(guān)系密切，如T 細胞耐受誘導的調(diào)節(jié)、CD4+的正向調(diào)節(jié)、T 細胞活化與分化的正向調(diào)節(jié)等生物過程，表明IDD 與T 細胞的調(diào)節(jié)有關(guān)，T 細胞參與IDD 的病理過程。而結(jié)果中的其他相關(guān)通路，如Fox O 信號通路、MAPK 信號通路與炎性反應及氧化應激關(guān)系密切［17-18］，鑒于炎性反應是IDD 或椎間盤再生機制的中心，其平衡是維持組織穩(wěn)態(tài)的關(guān)鍵［19］，抑制炎性反應可能有助于治療IDD［20］。

本研究使用免疫細胞浸潤分析的方法探索與IDD 具有強相關(guān)性的免疫細胞，結(jié)果發(fā)現(xiàn)，IDD患者與對照組的CD4 naive 表達差異性最為顯著。Naive T 細胞也稱非致敏T 細胞，在暴露于抗原刺激之前處于相對靜止的狀態(tài)。CD4 T 細胞是人體免疫系統(tǒng)中重要的免疫細胞［21］。最近的研究［22］表明，CD4 naive 的功能比以前認為的要多，并且它們在表型、分化階段、持久性、功能和解剖位置方面是多種多樣的。動物實驗［19］表明，椎間盤內(nèi)促炎治療可使血液和淋巴結(jié)中的T 細胞亞群增加。另一項研究［23］表明，在髓核突出的椎間盤中存在CD4 和CD8 T 細胞。在本研究中，CD4 naive 和CD4 T 細胞在免疫細胞浸潤分析后均存在顯著差異。Pearson 相關(guān)分析發(fā)現(xiàn)，UHMK1、ZFP36L2、ZCCHC3、ZBTB20和LINC00641 可能是IDD 免疫相關(guān)的關(guān)鍵靶點，即IDD 的特征基因。本研究還篩選了1 個與CD4 naive相關(guān)的ceRNA 子網(wǎng)絡，并將LINC00641 識別為關(guān)鍵因素，可將其作為IDD的潛在治療靶點。GSEA結(jié)果顯示，轉(zhuǎn)化生長因子信號通路可能是IDD 潛在的信號通路。上述關(guān)鍵基因的表達在GSE124272 數(shù)據(jù)集中得到了驗證。

目前大多數(shù)與LINC00641相關(guān)的研究都與惡性腫瘤有關(guān)［24］，也與糖尿病足的一些生物學過程有關(guān)［25］。在IDD 的研究中，Wang 等［26］發(fā)現(xiàn)LINC00641 可以直接與miR-153-3p 結(jié)合，作為miR-153-3p 的天然海綿，通過miRNA介導的沉默途徑靶向自噬相關(guān)基因5，并調(diào)節(jié)髓核細胞的自噬和死亡。本研究結(jié)果表明，LINC00641 可能通過免疫相關(guān)競爭性調(diào)節(jié)基因機制影響IDD 的發(fā)生與發(fā)展。

本研究發(fā)現(xiàn)的與CD4 naive 相關(guān)的mRNA 目前研究尚少，這些mRNA與腫瘤和免疫密切相關(guān)［27-29］。UHMK1 可能參與RNA 加工，也可能參與白細胞的形成［30］。ZFP36L2 在CD4 naive 中高表達，可抑制調(diào)節(jié)性T 細胞的功能［31］。本研究的免疫浸潤分析結(jié)果顯示，CD4 和CD4 naive 與IDD 密切相關(guān)，間接驗證了ZFP36L2 是IDD 的重要靶點之一。同時，ZBTB20 可抑制CD8+T 細胞免疫代謝［32］，并能促進巨噬細胞炎性反應，啟動先天免疫反應［33］。這些結(jié)果為IDD 的診斷和靶向治療提供了新思路。

本研究通過GSEA 來識別LINC00641 的關(guān)鍵調(diào)控途徑，發(fā)現(xiàn)與免疫相關(guān)的關(guān)鍵途徑是轉(zhuǎn)化生長因子信號通路、胞質(zhì)DNA 感應通路、溶酶體通路等。胞質(zhì)DNA 感應通路與先天免疫反應密切相關(guān)［34］，溶酶體是真核細胞的主要分解代謝區(qū)［35］，許多炎性疾病的根本原因可能是自噬-溶酶體途徑的功能障礙［36］。最近有研究［37］表明，溶酶體與椎間盤細胞的死亡有關(guān)。近年來，有研究［38］發(fā)現(xiàn)轉(zhuǎn)化生長因子在調(diào)節(jié)細胞生長分化和免疫功能方面具有重要作用。轉(zhuǎn)化生長因子信號通路對椎間盤的發(fā)育和生長至關(guān)重要［39］，但過度激活轉(zhuǎn)化生長因子信號對椎間盤是有害的，抑制異常的轉(zhuǎn)化生長因子信號可以延緩椎間盤的退行性變［40］。因此，通過lncRNA 抑制轉(zhuǎn)化生長因子信號通路來減少炎性反應可能是治療IDD 的有效方法。

本研究尚有一定的局限性，受現(xiàn)有數(shù)據(jù)庫中數(shù)據(jù)的限制，本研究選取的樣本量有限，以后將增加樣本量進一步研究探索；此外，本研究組未開展細胞實驗或動物實驗，未來將通過基礎(chǔ)實驗研究進一步驗證關(guān)鍵lncRNA 和mRNA 在IDD 中的作用。

總之，本研究在WGCNA 的基礎(chǔ)上，將ceRNA、IDD 和免疫調(diào)節(jié)聯(lián)系起來，構(gòu)建了一個與免疫調(diào)節(jié)相關(guān)的IDD ceRNA 網(wǎng)絡。在這個網(wǎng)絡中，lncRNA LINC00641 是一個關(guān)鍵靶點，它干預了mRNA（UHMK1、ZFP36L2、ZCCHC3、ZBTB20）的表達。同時，CD4 naive 表達減少，這可能與IDD 的發(fā)生機制密切相關(guān)。lncRNA LINC00641 可能通過轉(zhuǎn)化生長因子信號通路調(diào)節(jié)IDD。這些結(jié)果將有助于從ceRNA 的角度更好地探索免疫調(diào)節(jié)在IDD 發(fā)展中的作用機制，并為治療靶點的開發(fā)提供思路。