基于生物信息學的方法對骨肉瘤潛在樞紐基因的篩選及驗證

徐 磊,馬 濤

(皖南醫(yī)學院第一附屬醫(yī)院弋磯山醫(yī)院手足踝外科,安徽 蕪湖 241001)

骨肉瘤是最常見的原發(fā)性骨骼惡性腫瘤[1],占兒童惡性骨腫瘤的50%以上,是青少年癌癥死亡的第二大原因[2]。過去常用的治療方法是手術切除,但這種方法并沒有提高轉(zhuǎn)移性骨肉瘤患者的生存率。強化化療等現(xiàn)代聯(lián)合治療現(xiàn)已使骨肉瘤患者的五年生存率提高了35%～50%,然而,其高度的侵襲性和高轉(zhuǎn)移傾向使骨肉瘤的五年生存率仍然很低[3]。因此研究骨肉瘤發(fā)生發(fā)展機制探索治療新靶點顯得尤為重要。本研究基于GEO數(shù)據(jù)庫中的骨肉瘤相關數(shù)據(jù)集篩選及驗證骨肉瘤潛在樞紐基因,并通過TARGET數(shù)據(jù)庫中的骨肉瘤臨床信息對篩選出來的樞基因進行生存分析,以期為骨肉瘤的發(fā)病機制及發(fā)掘新的治療靶點提供思路。

1 資料與方法

1.1 基因表達譜芯片數(shù)據(jù)的獲取從GEO數(shù)據(jù) (http://www.ncbi.nlm.nih.gov/geo/)獲取與骨肉瘤相關的GSE16088及GSE36001數(shù)據(jù)集,分別基于GPL96平臺及GPL6102平臺。GSE16088數(shù)據(jù)集包含有骨肉瘤樣本14例,正常對照樣本9例。GSE36001數(shù)據(jù)集包含有骨肉瘤樣本19例,正常對照樣本6例。

1.2 數(shù)據(jù)預處理及DEGs篩選使用GEOquery包從GEO數(shù)據(jù)庫中下載 GSE16088 及GSE36001,并進行數(shù)據(jù)的預處理,然后通過箱式圖查看樣本標準化的情況,通過 PCA圖 以及 UMAP圖查看樣本分組間聚類情況,利用limma包進行兩組樣本的差異分析,以P<0.05、|logFC| ≥1.0為標準篩選骨肉瘤樣本及正常對照樣本的差異基因(differentially expressed genes,DEGS),使用ggplot2包繪制DEGs火山圖,用ComplexHeatmap包繪制熱圖。用韋恩圖(https://bioinformatics.psb.ugent.be/webtools/Venn/)獲取以上兩個數(shù)據(jù)集的共表達DEGs。

1.3 功能富集分析通常我們使用DAVID(https://david.ncifcrf.gov/)數(shù)據(jù)庫進行基因的富集分析[4]。打開DAVID在線工具,將GSE16088及GSE36001數(shù)據(jù)集共表達的155個DEGs上傳至列表后,得到GO分析包括生物學過程 (biological process,BP)、分子功能(molecular function,MF)、細胞成分(cellular component,CC)及 KEGG 通路的富集分析結(jié)果,P<0.05為差異具有統(tǒng)計學意義,使用R語言ggplot包繪制氣泡圖。

1.4 網(wǎng)絡構(gòu)建及樞紐基因的篩選PPI網(wǎng)絡可以幫助我們從相互作用層面識別參與骨肉瘤發(fā)展的樞紐基因和核心基因模塊。共表達DEGs的PPI信息是從String數(shù)據(jù)庫 (http://www.string-db.org/) 中獲得的(設置combined score>0.4)。我們將得到的PPI信息導入 Cytoscape軟件中進行可視化,并通過MCODE插件,以 degree cutoff=2、node score cutoff=0.2、k-core=2和max.depth=100為標準,篩選出 PPI 中的功能模塊。最后我們應用 Cytohubba插件,按Degree算法計算出PPI網(wǎng)絡中評分前6位的樞紐基因。

1.5 樞紐基因的驗證選擇GEO數(shù)據(jù)庫中的骨肉瘤相關數(shù)據(jù)集GSE33382,該數(shù)據(jù)集基于GPL10295平臺,包含有84個骨肉瘤樣本及3個成骨細胞對照樣本,驗證篩選出的6個樞紐基因在骨肉瘤樣本與對照樣本之間的差異表達。最后我們通TARGET數(shù)據(jù)庫骨肉瘤患者臨床隨訪資料及RNAseq相關數(shù)據(jù)對6個樞紐基因進行生存分析,評估樞紐基因?qū)侨饬龌颊叩念A后價值。P<0.05為差異有統(tǒng)計學意義。

2 結(jié)果

2.1 DEGs的篩選我們從GSE16088數(shù)據(jù)集共獲得了骨肉瘤樣本及正常對照樣本的2 649個DEGs,其中上調(diào)基因1 931個,下調(diào)基因718個。我們從GSE36001數(shù)據(jù)集共獲得了骨肉瘤樣本及正常對照樣本的754個DEGs,其中上調(diào)基因502個,下調(diào)基因252個。使用火山圖及熱圖可視化基因表達情況(見圖1、圖2)。韋恩圖獲取兩數(shù)據(jù)集共表達基因155個(見圖3)。

圖1 DEGs表達火山圖

圖2 DEGs高低表達各top20基因熱圖

圖3 GSE16088與GSE36001 韋恩圖

2.2 DEGs的GO富集分析和 KEGG信號通路分析GO富集分析的結(jié)果因GO分類和DEGs的表達變化而不同。對GSE16088及GSE36001數(shù)據(jù)集共表達的155個DEGs進行GO和KEGG分析。GO分析顯示共表達DEGs生物學功能主要參與調(diào)控中性粒細胞脫顆粒、細胞外結(jié)構(gòu)組織、腫瘤壞死因子產(chǎn)生的正向調(diào)節(jié)、骨吸收、含有膠原的細胞外基質(zhì)、溶酶體腔、膜區(qū)、分泌顆粒膜等。KEGG分析顯示共表達DEGs主要和 溶酶體、破骨細胞分化、利什曼病、細胞粘附分子、Th1及Th2細胞分化等有關(見表1、圖4)。

表1 共表達 DEGs GO富集分析及KEGG信號通路分析

圖4 GSE16088及GSE36001共表達DEGs GO+KEGG富集分析柱狀圖

2.3 DEGs的PPI網(wǎng)絡分析及樞紐基因的篩選通過String在線工具我們獲得了combined score>0.4的PPI信息及網(wǎng)絡圖,將PPI信息導入Cytoscape 軟件進行可視化獲得127節(jié)點396邊的PPI網(wǎng)絡圖(見圖5),利用Cytoscape內(nèi)的MCODE插件篩選出了得分為8.44,由10個節(jié)點38條邊組成的功能模塊(見圖6),利用cytohubba插件我們篩選出了前6位的樞紐基因,樞紐基因分別是酪氨酸激酶結(jié)合蛋白(tyrosine kinase binding protein,TYROBP)、基質(zhì)金屬蛋白酶9(matrix metalloproteinases9,MMP9)、集落刺激因子1受體(colony-stimulating factor 1 receptor,CSF1R)、細胞色素B-245 β鏈(cytochrome B-245 beta chain,CYBB)、分泌型磷蛋白1(secreted phosphoprotein 1,SPP1)、血小板內(nèi)皮細胞粘附分子1(platelet endothelial cell adhesion molecule-1,PECAM1)(見圖7)。這些樞紐基因在骨肉瘤樣本中均表達上調(diào)。

圖5 GSE16088及GSE36001共表達共表達 DEGs的PPI網(wǎng)絡

圖6 PPI網(wǎng)絡中篩選出的功能模塊

圖7 從功能模塊中篩選前6位樞紐基因

2.4 樞紐基因的驗證及與骨肉瘤患者預后分析選擇GEO數(shù)據(jù)庫中與骨肉瘤相關的數(shù)據(jù)集GSE33382,分析骨肉瘤樣本與正常對照樣本之間TYROBP、MMP9、CSF1R、CYBB、SPP1、PECAM1的差異表達水平,這6個樞紐基因在骨肉瘤樣本中均顯著上調(diào)表達(P<0.05,見圖8)。最后我們通TARGET數(shù)據(jù)庫骨肉瘤患者臨床隨訪資料及RNAseq相關數(shù)據(jù)對6個樞紐基因進行生存分析,發(fā)現(xiàn)TYROBP、CSF1R、CYBB高表達與骨肉瘤的總生存率成正相關(P<0.05),余下三個基因差異無統(tǒng)計學意義(見圖9)。

圖8 樞紐基因在GSE33382數(shù)據(jù)集中骨肉瘤及正常對照樣本的表達驗證

圖9 TYROBP、CYBB、CSF1R高低表達生存曲線

3 討論

本研究通過分析GEO數(shù)據(jù)庫中與骨肉瘤相關的數(shù)據(jù)集GSE16088及GSE36001得到共表達DEGs 155 個。GO分析顯示共表達DEGs生物學功能主要參與調(diào)控中性粒細胞脫顆粒、細胞外結(jié)構(gòu)組織、腫瘤壞死因子產(chǎn)生的正向調(diào)節(jié)、骨吸收、含有膠原的細胞外基質(zhì)、溶酶體腔、膜區(qū)、分泌顆粒膜等,KEGG分析顯示共表達DEGs主要和 溶酶體、破骨細胞分化、利什曼病、細胞粘附分子、Th1及Th2細胞分化等有關。應用STRING數(shù)據(jù)庫和Cytoscape構(gòu)建了共表達DEGs的PPI網(wǎng)絡圖,并篩選出6個樞紐基因,分別是 TYROBP、MMP9、CSF1R、CYBB、SPP1、PECAM1。經(jīng)GSE33382數(shù)據(jù)集驗證,6個樞紐基因在骨肉瘤樣本中均高表達。且經(jīng)過生存分析顯示TYROBP、CSF1R、CYBB高表達與骨肉瘤患者的有利預后相關。

TYROBP是一種基于免疫受體酪氨酸的激活基序,主要在自然殺傷細胞和骨髓細胞中表達,可與幾種免疫受體結(jié)合激活多種信號通路[5]。Liang等[6]研究發(fā)現(xiàn)TYROBP主要富集于自然殺傷細胞介導的細胞毒性和破骨細胞分化,可導致腫瘤細胞凋亡并促進破骨細胞分化,從而引起腫瘤周圍的骨吸收,TYROBP的低表達可促進骨肉瘤的發(fā)生和進展。我們研究發(fā)現(xiàn)TYROBP的高表達與骨肉瘤的有利預后相關,與此研究的結(jié)果是一致的。MMP9又稱為IV型膠原酶及明膠酶B,由 Wilhelm于1989年發(fā)現(xiàn),是最大的基質(zhì)金屬蛋白酶MMP[7]。目前已有很多研究證實了MMP9與骨肉瘤的增殖、侵襲及轉(zhuǎn)移等有關[8-9]。CSF1R屬于III型蛋白酪氨酸激酶受體家族,與CSF1或最近發(fā)現(xiàn)的配體IL-34的結(jié)合可誘導受體同源二聚化并隨后激活受體信號傳導[10]。Wen等研究發(fā)現(xiàn)骨肉瘤細胞系及骨肉瘤活檢標本中有CSF1R的表達,并且高表達CSF1R有促進骨肉瘤細胞增殖、存活及轉(zhuǎn)移的作用,其機制與巨噬細胞表達的CSF-1R的促腫瘤發(fā)生作用相似[11]。既往研究表明CSF1R+巨噬細胞(M2型)可導致免疫抑制從而有促進腫瘤發(fā)生的作用[12]。腫瘤微環(huán)境在肉瘤形成和免疫治療中起著重要作用,因此,越來越多的研究集中于骨肉瘤的腫瘤微環(huán)境。Song等[13]研究發(fā)現(xiàn)CD4+Th1細胞和CD68+巨噬細胞的存在可能預示著骨肉瘤的有利預后,而CSF1R 與 CD4/CD68之間存在正相關性,間接顯示了CSFIR與骨肉瘤的有利預后之間的聯(lián)系,但具體的機制目前還不清楚。CYBB也稱為CYBB/NOX2,為NADPH氧化酶蛋白家族中的一員。它是第一個在吞噬細胞內(nèi)體中被鑒定并高表達的NADPH氧化酶,已知具有多種功能,對抗菌宿主防御至關重要[14]。目前已證實CYBB在包括肺癌、肝癌在內(nèi)的部分人類腫瘤中發(fā)揮重要作用[15-16]。CYBB的mRNA水平在骨肉瘤細胞系中過表達,抑制CYBB表達可誘導骨肉瘤細胞凋亡[17]。Lin等[18]發(fā)現(xiàn)與非腫瘤性骨組織相比,CYBB 在骨肉瘤中顯著過表達,并且與骨肉瘤的不良預后相關。SPP1,又稱為骨橋蛋白,是一種整合素結(jié)合磷酸化糖蛋白,據(jù)報道在多種腫瘤細胞中差異表達,和腫瘤的增殖、侵襲、遷移、血管生成有關[19]。LAMP3可能通過調(diào)節(jié)SPP1下游的信號傳導參與骨肉瘤的侵襲和轉(zhuǎn)移[20]。Dalla-Torre等[21]研究發(fā)現(xiàn)SPP1高表達的骨肉瘤患者預后好。PECAM1在內(nèi)皮細胞緊密連接、血小板和免疫細胞大量表達,已證實抑制 PECAM1 可減少多種人類疾病的炎癥反應[22]。目前有關PECAM1與骨肉瘤發(fā)病機制的研究較少。邢江浩等[23]發(fā)現(xiàn)骨肉瘤細胞分泌的外泌體可進入血管內(nèi)皮細胞,導致血管生長能力及PECAM1 表達上調(diào)。鑒于PECAM1是血管內(nèi)皮細胞較為特征性的標志物,可能在骨肉瘤的侵襲及轉(zhuǎn)移中起作用。

綜上所述,本研究通過對GEO數(shù)據(jù)庫中GSE16088及GSE36001的生物信息學分析,篩選出了6個與骨肉瘤 相關的樞紐基因TYROBP、MMP9、CSF1R、CYBB、SPP1、PECAM1,其中TYROBP、CSF1R、CYBB高表達與骨肉瘤患者的有利預后相關,為骨肉瘤的發(fā)病的分子機制研究、預后判斷及治療提供了參考。