基于生物信息学的方法对骨肉瘤潜在枢纽基因的筛选及验证

徐 磊,马 涛

(皖南医学院第一附属医院弋矶山医院手足踝外科,安徽 芜湖 241001)

骨肉瘤是最常见的原发性骨骼恶性肿瘤[1],占儿童恶性骨肿瘤的50%以上,是青少年癌症死亡的第二大原因[2]。过去常用的治疗方法是手术切除,但这种方法并没有提高转移性骨肉瘤患者的生存率。强化化疗等现代联合治疗现已使骨肉瘤患者的五年生存率提高了35%～50%,然而,其高度的侵袭性和高转移倾向使骨肉瘤的五年生存率仍然很低[3]。因此研究骨肉瘤发生发展机制探索治疗新靶点显得尤为重要。本研究基于GEO数据库中的骨肉瘤相关数据集筛选及验证骨肉瘤潜在枢纽基因,并通过TARGET数据库中的骨肉瘤临床信息对筛选出来的枢基因进行生存分析,以期为骨肉瘤的发病机制及发掘新的治疗靶点提供思路。

1 资料与方法

1.1 基因表达谱芯片数据的获取从GEO数据 (http://www.ncbi.nlm.nih.gov/geo/)获取与骨肉瘤相关的GSE16088及GSE36001数据集,分别基于GPL96平台及GPL6102平台。GSE16088数据集包含有骨肉瘤样本14例,正常对照样本9例。GSE36001数据集包含有骨肉瘤样本19例,正常对照样本6例。

1.2 数据预处理及DEGs筛选使用GEOquery包从GEO数据库中下载 GSE16088 及GSE36001,并进行数据的预处理,然后通过箱式图查看样本标准化的情况,通过 PCA图 以及 UMAP图查看样本分组间聚类情况,利用limma包进行两组样本的差异分析,以P<0.05、|logFC| ≥1.0为标准筛选骨肉瘤样本及正常对照样本的差异基因(differentially expressed genes,DEGS),使用ggplot2包绘制DEGs火山图,用ComplexHeatmap包绘制热图。用韦恩图(https://bioinformatics.psb.ugent.be/webtools/Venn/)获取以上两个数据集的共表达DEGs。

1.3 功能富集分析通常我们使用DAVID(https://david.ncifcrf.gov/)数据库进行基因的富集分析[4]。打开DAVID在线工具,将GSE16088及GSE36001数据集共表达的155个DEGs上传至列表后,得到GO分析包括生物学过程 (biological process,BP)、分子功能(molecular function,MF)、细胞成分(cellular component,CC)及 KEGG 通路的富集分析结果,P<0.05为差异具有统计学意义,使用R语言ggplot包绘制气泡图。

1.4 网络构建及枢纽基因的筛选PPI网络可以帮助我们从相互作用层面识别参与骨肉瘤发展的枢纽基因和核心基因模块。共表达DEGs的PPI信息是从String数据库 (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网络中评分前6位的枢纽基因。

1.5 枢纽基因的验证选择GEO数据库中的骨肉瘤相关数据集GSE33382,该数据集基于GPL10295平台,包含有84个骨肉瘤样本及3个成骨细胞对照样本,验证筛选出的6个枢纽基因在骨肉瘤样本与对照样本之间的差异表达。最后我们通TARGET数据库骨肉瘤患者临床随访资料及RNAseq相关数据对6个枢纽基因进行生存分析,评估枢纽基因对骨肉瘤患者的预后价值。P<0.05为差异有统计学意义。

2 结果

2.1 DEGs的筛选我们从GSE16088数据集共获得了骨肉瘤样本及正常对照样本的2 649个DEGs,其中上调基因1 931个,下调基因718个。我们从GSE36001数据集共获得了骨肉瘤样本及正常对照样本的754个DEGs,其中上调基因502个,下调基因252个。使用火山图及热图可视化基因表达情况(见图1、图2)。韦恩图获取两数据集共表达基因155个(见图3)。

图1 DEGs表达火山图

图2 DEGs高低表达各top20基因热图

图3 GSE16088与GSE36001 韦恩图

2.2 DEGs的GO富集分析和 KEGG信号通路分析GO富集分析的结果因GO分类和DEGs的表达变化而不同。对GSE16088及GSE36001数据集共表达的155个DEGs进行GO和KEGG分析。GO分析显示共表达DEGs生物学功能主要参与调控中性粒细胞脱颗粒、细胞外结构组织、肿瘤坏死因子产生的正向调节、骨吸收、含有胶原的细胞外基质、溶酶体腔、膜区、分泌颗粒膜等。KEGG分析显示共表达DEGs主要和 溶酶体、破骨细胞分化、利什曼病、细胞粘附分子、Th1及Th2细胞分化等有关(见表1、图4)。

表1 共表达 DEGs GO富集分析及KEGG信号通路分析

图4 GSE16088及GSE36001共表达DEGs GO+KEGG富集分析柱状图

2.3 DEGs的PPI网络分析及枢纽基因的筛选通过String在线工具我们获得了combined score>0.4的PPI信息及网络图,将PPI信息导入Cytoscape 软件进行可视化获得127节点396边的PPI网络图(见图5),利用Cytoscape内的MCODE插件筛选出了得分为8.44,由10个节点38条边组成的功能模块(见图6),利用cytohubba插件我们筛选出了前6位的枢纽基因,枢纽基因分别是酪氨酸激酶结合蛋白(tyrosine kinase binding protein,TYROBP)、基质金属蛋白酶9(matrix metalloproteinases9,MMP9)、集落刺激因子1受体(colony-stimulating factor 1 receptor,CSF1R)、细胞色素B-245 β链(cytochrome B-245 beta chain,CYBB)、分泌型磷蛋白1(secreted phosphoprotein 1,SPP1)、血小板内皮细胞粘附分子1(platelet endothelial cell adhesion molecule-1,PECAM1)(见图7)。这些枢纽基因在骨肉瘤样本中均表达上调。

图5 GSE16088及GSE36001共表达共表达 DEGs的PPI网络

图6 PPI网络中筛选出的功能模块

图7 从功能模块中筛选前6位枢纽基因

2.4 枢纽基因的验证及与骨肉瘤患者预后分析选择GEO数据库中与骨肉瘤相关的数据集GSE33382,分析骨肉瘤样本与正常对照样本之间TYROBP、MMP9、CSF1R、CYBB、SPP1、PECAM1的差异表达水平,这6个枢纽基因在骨肉瘤样本中均显著上调表达(P<0.05,见图8)。最后我们通TARGET数据库骨肉瘤患者临床随访资料及RNAseq相关数据对6个枢纽基因进行生存分析,发现TYROBP、CSF1R、CYBB高表达与骨肉瘤的总生存率成正相关(P<0.05),余下三个基因差异无统计学意义(见图9)。

图8 枢纽基因在GSE33382数据集中骨肉瘤及正常对照样本的表达验证

图9 TYROBP、CYBB、CSF1R高低表达生存曲线

3 讨论

本研究通过分析GEO数据库中与骨肉瘤相关的数据集GSE16088及GSE36001得到共表达DEGs 155 个。GO分析显示共表达DEGs生物学功能主要参与调控中性粒细胞脱颗粒、细胞外结构组织、肿瘤坏死因子产生的正向调节、骨吸收、含有胶原的细胞外基质、溶酶体腔、膜区、分泌颗粒膜等,KEGG分析显示共表达DEGs主要和 溶酶体、破骨细胞分化、利什曼病、细胞粘附分子、Th1及Th2细胞分化等有关。应用STRING数据库和Cytoscape构建了共表达DEGs的PPI网络图,并筛选出6个枢纽基因,分别是 TYROBP、MMP9、CSF1R、CYBB、SPP1、PECAM1。经GSE33382数据集验证,6个枢纽基因在骨肉瘤样本中均高表达。且经过生存分析显示TYROBP、CSF1R、CYBB高表达与骨肉瘤患者的有利预后相关。

TYROBP是一种基于免疫受体酪氨酸的激活基序,主要在自然杀伤细胞和骨髓细胞中表达,可与几种免疫受体结合激活多种信号通路[5]。Liang等[6]研究发现TYROBP主要富集于自然杀伤细胞介导的细胞毒性和破骨细胞分化,可导致肿瘤细胞凋亡并促进破骨细胞分化,从而引起肿瘤周围的骨吸收,TYROBP的低表达可促进骨肉瘤的发生和进展。我们研究发现TYROBP的高表达与骨肉瘤的有利预后相关,与此研究的结果是一致的。MMP9又称为IV型胶原酶及明胶酶B,由 Wilhelm于1989年发现,是最大的基质金属蛋白酶MMP[7]。目前已有很多研究证实了MMP9与骨肉瘤的增殖、侵袭及转移等有关[8-9]。CSF1R属于III型蛋白酪氨酸激酶受体家族,与CSF1或最近发现的配体IL-34的结合可诱导受体同源二聚化并随后激活受体信号传导[10]。Wen等研究发现骨肉瘤细胞系及骨肉瘤活检标本中有CSF1R的表达,并且高表达CSF1R有促进骨肉瘤细胞增殖、存活及转移的作用,其机制与巨噬细胞表达的CSF-1R的促肿瘤发生作用相似[11]。既往研究表明CSF1R+巨噬细胞(M2型)可导致免疫抑制从而有促进肿瘤发生的作用[12]。肿瘤微环境在肉瘤形成和免疫治疗中起着重要作用,因此,越来越多的研究集中于骨肉瘤的肿瘤微环境。Song等[13]研究发现CD4+Th1细胞和CD68+巨噬细胞的存在可能预示着骨肉瘤的有利预后,而CSF1R 与 CD4/CD68之间存在正相关性,间接显示了CSFIR与骨肉瘤的有利预后之间的联系,但具体的机制目前还不清楚。CYBB也称为CYBB/NOX2,为NADPH氧化酶蛋白家族中的一员。它是第一个在吞噬细胞内体中被鉴定并高表达的NADPH氧化酶,已知具有多种功能,对抗菌宿主防御至关重要[14]。目前已证实CYBB在包括肺癌、肝癌在内的部分人类肿瘤中发挥重要作用[15-16]。CYBB的mRNA水平在骨肉瘤细胞系中过表达,抑制CYBB表达可诱导骨肉瘤细胞凋亡[17]。Lin等[18]发现与非肿瘤性骨组织相比,CYBB 在骨肉瘤中显著过表达,并且与骨肉瘤的不良预后相关。SPP1,又称为骨桥蛋白,是一种整合素结合磷酸化糖蛋白,据报道在多种肿瘤细胞中差异表达,和肿瘤的增殖、侵袭、迁移、血管生成有关[19]。LAMP3可能通过调节SPP1下游的信号传导参与骨肉瘤的侵袭和转移[20]。Dalla-Torre等[21]研究发现SPP1高表达的骨肉瘤患者预后好。PECAM1在内皮细胞紧密连接、血小板和免疫细胞大量表达,已证实抑制 PECAM1 可减少多种人类疾病的炎症反应[22]。目前有关PECAM1与骨肉瘤发病机制的研究较少。邢江浩等[23]发现骨肉瘤细胞分泌的外泌体可进入血管内皮细胞,导致血管生长能力及PECAM1 表达上调。鉴于PECAM1是血管内皮细胞较为特征性的标志物,可能在骨肉瘤的侵袭及转移中起作用。

综上所述,本研究通过对GEO数据库中GSE16088及GSE36001的生物信息学分析,筛选出了6个与骨肉瘤 相关的枢纽基因TYROBP、MMP9、CSF1R、CYBB、SPP1、PECAM1,其中TYROBP、CSF1R、CYBB高表达与骨肉瘤患者的有利预后相关,为骨肉瘤的发病的分子机制研究、预后判断及治疗提供了参考。