文章信息
- 裴忠霞, 孙晶
- PEI Zhongxia, SUN Jing
- 基于GEO和TARGET数据库分析骨肉瘤转移和预后的关键基因
- Analysis of hub genes for predicting metastasis and prognosis of osteosarcoma in GEO and TARGET databases
- 中国医科大学学报, 2023, 52(11): 1018-1024
- Journal of China Medical University, 2023, 52(11): 1018-1024
-
文章历史
- 收稿日期:2022-12-13
- 网络出版时间:2023-11-02 11:35:03
骨肉瘤常见于儿童、青少年以及60岁以上人群,具有预后差、致残率高的特点[1]。骨肉瘤具有较高的局部侵袭和早期转移倾向[2]。目前,有关骨肉瘤患者的治疗方法主要包括手术、放疗、化疗和新辅助化疗等,很大程度上改善了骨肉瘤患者的预后[3],但转移性和复发性骨肉瘤的预后仍不理想,5年生存率仅为50%[4]。因此,探索骨肉瘤转移的发病机制,寻找相关预后指标,有助于为转移性骨肉瘤患者的治疗提供新的靶点。
近年来,随着高通量测序技术的发展,基于生物信息学方法进行的大数据挖掘分析被广泛用于探索肿瘤的诊断和预后相关的生物标志物[5]。本研究使用基因表达综合数据库(Gene Expression Omnibus,GEO)和TARGET等公共数据库筛选和分析骨肉瘤转移相关的差异表达基因(differentially expressed genes,DEGs),探讨其对骨肉瘤患者预后的影响,为骨肉瘤转移的分子诊断研究和临床治疗提供新的理论依据。
1 材料与方法 1.1 数据收集从GEO(https://www.ncbi.nlm.nih.gov/geo/)中筛选骨肉瘤转移相关的数据集,获取符合研究的基因芯片数据集GSE21257,芯片平台为GPL10295。该数据集包含5年内未发生转移的骨肉瘤患者样本19例和5年内发生转移的骨肉瘤患者样本34例,根据转移情况将样本分为未转移组和转移组。
1.2 筛选2组样本之间的DEGs采用GEO2R(https://www.ncbi.nlm.nih.gov/geo/geo2r/)在线分析工具检测未转移组与转移组样本之间的DEGs。具体筛选标准为调整后P < 0.05且log2|FC| > 1。DEGs结果通过火山图和热图分别表示。
1.3 DEGs的相关功能富集分析通过DAVID数据库(https://david.ncifcrf.gov/)对DEGs进行基因本体(GeneOntology,GO)功能注释和京都基因与基因组数据库(Kyoto Encyclopedia of Genesand Genomes,KEGG)分析。其中GO分析包括生物过程(biologicalprocess,BP)、细胞成分(cellularcomponent,CC)和分子功能(molecularfunction,MF)。以调整后P < 0.05为条件,筛选有统计学差异的蛋白质功能和相关信号通路。
1.4 构建蛋白质-蛋白质相互作用(protein-proteininteration,PPI)网络及关键基因的筛选为了明确DEGs的蛋白质之间相互作用关系,采用STRING(https://string-db.org)数据库建立DEGs的PPI网络。筛选最低交互分数 > 0.4的相互作用蛋白,随后将得到的蛋白质数据导入到Cytoscape软件(v3.9.1),按照Betweenness Centrality(BC)数值的大小绘制PPI同心圆,随后使用cytoHubba插件按照度值算法选取PPI图中排名前5的关键基因用于骨肉瘤的生存等相关研究。
1.5 采用TARGET数据库分析骨肉瘤患者的生存差异TARGET(https://ocg.cancer.gov/programs/target/)是目前最大的儿童肿瘤数据库,该数据库纳入了87例骨肉瘤患者完整的生存数据(其中1例数据不完整,下载时间为2022年11月10日)。从TARGET数据库中获得了骨肉瘤患者的基因表达和临床数据矩阵。多基因相关性图通过R软件包pheatmap进行展示,基因的单因素Cox分析通过R软件包forestplot进行实现。将关键基因的表达按照高低进行分组,通过R软件包survival分析关键基因对骨肉瘤患者预后的影响。
1.6 统计学分析应用GraphPad Prism 9.4.1对数据进行统计分析和绘图。多组间的比较采用单因素方差分析,组间两两比较采用LSD-t检验,采用单因素Cox回归分析关键基因对患者预后的影响,采用Spearman相关性分析检验关键基因之间的相关性,log-rank用于检验Kaplan-Meier(K-M)曲线生存分析,比较2组之间的生存差异。P < 0.05为差异有统计学意义。
2 结果 2.1 DEGs的分析和筛选GEO2R在线分析共得到55个DEGs。其中下调基因22个,上调基因33个。分别绘制DEGs的火山图(图 1A)和前10个上调差异基因与前10个下调差异基因的热图(图 1B)。
![]() |
A, volcano plot on which the red dots indicate significantly differentially up regulated genes, and the blue dots indicate significantly differentially down regulated genes. B, heat map. 图 1 GSE21257芯片数据中DEGs的火山图和热图 Fig.1 Volcano plot and heat map of DEGs in the GSE21257 gene microarray |
2.2 GO功能富集和KEGG通路富集分析
通过DAVID在线数据分析软件对DEGs进行GO和KEGG分析,按照调整后P值大小排序,分别选取排名前3位BP、CC、MF和KEGG通路,见表 1。功能注释结果显示,BP方面,差异基因主要参与抗原加工和外源肽抗原通过MHC Ⅱ类呈递,肽抗原与MHCⅡ类蛋白复合体组装,免疫球蛋白的产生参与免疫球蛋白介导的免疫反应;CC方面,差异基因主要参与MHC Ⅱ类蛋白复合体,膜泡运输,内质网膜腔侧的整体成分;MF方面,差异基因主要参与MHC Ⅱ类蛋白复合体结合,MHC Ⅱ类受体活性,肽抗原结合。KEGG通路富集分析发现,差异基因主要参与黄色金葡萄球菌感染、哮喘和Ⅰ型糖尿病。
Category | Term | Count | Adjusted P value |
GO_BP | Antigen processing and presentation of exogenous peptide antigen via MHC classⅡ | 11 | 1.28E-16 |
GO_BP | Peptide antigen assembly with MHC class Ⅱ protein complex | 9 | 5.55E-15 |
GO_BP | Immunoglobulin production involved in immunoglobulin mediated immune response | 9 | 6.97E-15 |
GO_CC | MHC class Ⅱ protein complex | 10 | 1.36E-15 |
GO_CC | Transport vesicle membrane | 8 | 9.24E-10 |
GO_CC | Integral component of lumenal side of endoplasmic reticulum membrane | 7 | 5.70E-09 |
GO_MF | MHC class Ⅱ protein complex binding | 10 | 4.07E-15 |
GO_MF | MHC class Ⅱ receptor activity | 7 | 4.95E-11 |
GO_MF | Peptide antigen binding | 5 | 1.75E-04 |
KEGG_pathway | Staphylococcus aureus infection | 14 | 9.73E-14 |
KEGG_pathway | Asthma | 10 | 1.21E-12 |
KEGG_pathway | TypeⅠdiabetes mellitus | 10 | 2.15E-11 |
2.3 PPI分析与关键基因筛选
将GSE21257数据集的55个DEGs导入到STRING数据库,设定筛选阈值为0.4,绘制PPI网络图,共得到53个节点和178条边,平均节点度为6.72,见图 2A。随后,将STRING数据库分析的蛋白互作网络关系文件结果导入Cytoscape软件进行进一步分析,并以BC值大小绘制同心圆。随后用cytoHubba插件,根据度值得出排名前5位的关键基因为TYROBP、ITGB2、C1QB、C1QC、CD74,5个关键基因在骨肉瘤转移样本中均下调,见图 2B。
![]() |
A, PPI network diagram; B, top five differentially expressed genes in the PPI network. 图 2 DEGs的PPI网络图和5个关键基因 Fig.2 PPI network diagram of differentially expressed genes and five key genes |
2.4 关键基因在骨肉瘤样本中的相关性分析和生存分析
将TARGET数据库骨肉瘤样本基因表达矩阵和临床信息数据下载后,通过R软件包pheatmap分析关键基因在骨肉瘤样本中的相关性,结果显示,5个关键基因TYROBP、ITGB2、C1QB、C1QC、CD74在骨肉瘤样本中均呈正相关(图 3A)。通过分析TARGET数据库骨肉瘤患者的生存模型,发现关键基因TYROBP、C1QB和CD74的表达与骨肉瘤患者的预后有统计学差异,见图 3B(P < 0.05),此外,CD74与骨肉瘤患者的无病生存期(disease free survival,DFS)相关,见图 3C(P < 0.05)。进一步进行K-M生存曲线分析发现,高表达TYROBP、C1QB和CD74组总生存期(overall survival,OS)较长(P < 0.05),CD74基因的上调使骨肉瘤患者有良好的DFS(P < 0.05),预后较好,见图 4。
![]() |
A, correlation map of key genes; B, forest maps describing P value, risk coefficient hazard ratio for univariate Cox analysis of overall survival; C, forest maps describing P value, risk coefficient hazard ratio for univariate Cox analysis of disease-free survival. 图 3 关键基因在骨肉瘤样本中的相关性和基因单因素Cox分析 Fig.3 Correlation and gene univariate Cox analysis of key genes in osteosarcoma samples |
![]() |
A, K-M survival curve of OS distribution of TYROBP from TARGET database; B, K-M survival curve of OS distribution of CD74 from TARGET database; C, K-M survival curve of OS distribution of C1QB from TARGET database; D, K-M survival curve of DFS distribution of CD74 from TARGET database. 图 4 与预后相关的关键基因的生存曲线分析 Fig.4 Survival curve analysis of key genes associated with prognosis |
3 讨论
骨肉瘤起源于骨髓中的一个单细胞,最终会产生一个多克隆的、异质性的肿瘤块,被认为是分子畸变方面最复杂的癌症之一[6]。如果不及时治疗,90%的骨肉瘤患者将死于肺转移[7]。骨肉瘤最常见于股骨远端,其次为胫骨近端和肱骨近端[8]。最初针对骨肉瘤的治疗主要是外科手术,但致残率较高,严重影响患者的生活质量且早期极易发生肺部转移[9]。20世纪初,随着化疗技术的不断更新,手术联合化疗使骨肉瘤患者的5年生存率提高至52% [9]。然而,由于骨肉瘤的恶性程度高易转移的特点,患者的OS仍远不能令人满意,尤其是晚期骨肉瘤患者。因此,一些针对骨肉瘤发展过程的新疗法,如抗血管生成药物和免疫疗法已被应用于骨肉瘤的临床治疗,但其效性和潜在机制仍不清楚[4]。深入了解骨肉瘤发生发展相关的分子病理机制,有助于筛选骨肉瘤早期诊断、靶向及免疫治疗和预后分析的关键分子或生物标志物。
本研究基于生物信息学分析方法,首先将GSE21257芯片数据集[10] 5年内未发生转移的骨肉瘤患者与5年内发生转移的骨肉瘤患者的全基因组基因表达谱进行比较,筛选出55个DEGs。通过GO功能富集分析和KEGG信号通路分析,发现差异基因主要参与MHC Ⅱ类蛋白复合体组装,免疫球蛋白介导的免疫反应,MHC Ⅱ类受体活性,肽抗原结合等。这一结果提示参与抗原呈递基因的失调可能是骨肉瘤发生的早期事件,MHC Ⅱ只在抗原呈递细胞表面表达,如巨噬细胞、树突状细胞和B细胞[11]。抗原呈递细胞通过将MHC Ⅱ与多肽结合,向辅助T细胞呈现外源性多肽或内源性多肽,从而告知机体正在被入侵[12]。信号通路方面主要参与金黄色葡萄球菌感染、哮喘和Ⅰ型糖尿病等。与本研究结果一致,TUOHY等[13]发现金黄色葡萄球菌感染会上调骨肉瘤患者炎症免疫反应,从而抵消骨肉瘤引起的免疫抑制,为优化炎症刺激以触发抗骨肉瘤巨噬细胞反应提供了一种潜在的治疗策略。
利用PPI分析对差异基因进一步筛选得到5个骨肉瘤转移关键基因:TYROBP、ITGB2、C1QB、C1QC和CD74。研究[14-15]发现TYROBP和ITGB2的高表达与骨肉瘤患者更好的OS相关,是预测骨肉瘤的独立预后因素。另有研究[16-17]表明关键基因C1QB、C1QC和CD74可能是骨肉瘤免疫微环境重构的潜在指标,有助于预测骨肉瘤患者的预后。通过分析TARGET数据库中骨肉瘤患者的临床数据,发现关键基因在骨肉瘤样本中均呈正相关,其中TYROBP、C1QB和CD74的高表达与骨肉瘤患者良好的OS相关,此外,CD74的高表达提示骨肉瘤患者的DFS较长。
综上所述,本研究通过生物信息方法筛选出骨肉瘤转移和预后相关的关键基因TYROBP、C1QB和CD74,在5年内发生骨肉瘤转移的样本中均呈现低表达,可能与骨肉瘤的发生和发展相关。这些结果表明TYROBP、C1QB和CD74在骨肉瘤患者良好的预后中起到关键作用。但考虑到上述研究为数据挖掘形式,具体的作用机制尚需进一步的基础和临床实验研究证实。
[1] |
SIEGEL RL, MILLER KD, JEMAL A. Cancer statistics, 2018[J]. CA A Cancer J Clin, 2018, 68(1): 7-30. DOI:10.3322/caac.21442 |
[2] |
BERHE S, DANZER E, MEYERS P, et al. Unusual abdominal metastases in osteosarcoma[J]. J Pediatr Surg Case Rep, 2018, 28: 13-16. DOI:10.1016/j.epsc.2017.09.022 |
[3] |
BENJAMIN RS. Adjuvant and neoadjuvant chemotherapy for osteosarcoma: a historical perspective[J]. Adv Exp Med Biol, 2020, 1257: 1-10. DOI:10.1007/978-3-030-43032-0_1 |
[4] |
CHEN CL, XIE L, REN TT, et al. Immunotherapy for osteosarcoma: fundamental mechanism, rationale, and recent breakthroughs[J]. Cancer Lett, 2021, 500: 1-10. DOI:10.1016/j.canlet.2020.12.024 |
[5] |
WANG YH, ZHAO Y, BOLLAS A, et al. Nanopore sequencing technology, bioinformatics and applications[J]. Nat Biotechnol, 2021, 39(11): 1348-1365. DOI:10.1038/s41587-021-01108-x |
[6] |
BROWN HK, TELLEZ-GABRIEL M, HEYMANN D. Cancer stem cells in osteosarcoma[J]. Cancer Lett, 2017, 386: 189-195. DOI:10.1016/j.canlet.2016.11.019 |
[7] |
RITTER J, BIELACK SS. Osteosarcoma[J]. Ann Oncol, 2010, 21: vii320-vii325. DOI:10.1093/annonc/mdq276 |
[8] |
CHO WH, SONG WS, JEON DG, et al. Differential presentations, clinical courses, and survivals of osteosarcomas of the proximal humerus over other extremity locations[J]. Ann Surg Oncol, 2010, 17(3): 702-708. DOI:10.1245/s10434-009-0825-6 |
[9] |
BELAYNEH R, FOURMAN MS, BHOGAL S, et al. Update on osteosarcoma[J]. Curr Oncol Rep, 2021, 23(6): 71. DOI:10.1007/s11912-021-01053-7 |
[10] |
BUDDINGH EP, KUIJJER ML, DUIM RAJ, et al. Tumor-infiltrating macrophages are associated with metastasis suppression in high-grade osteosarcoma: a rationale for treatment with macrophage activating agents[J]. Clin Cancer Res, 2011, 17(8): 2110-2119. DOI:10.1158/1078-0432.CCR-10-2047 |
[11] |
MORI L, LEPORE M, DE LIBERO G. The immunology of CD1-and MR1-restricted T cells[J]. Annu Rev Immunol, 2016, 34: 479-510. DOI:10.1146/annurev-immunol-032414-112008 |
[12] |
LIU JW, WU SY, XIE XY, et al. Identification of potential crucial genes and key pathways in osteosarcoma[J]. Hereditas, 2020, 157(1): 1-13. DOI:10.1186/s41065-020-00142-0 |
[13] |
TUOHY JL, SOMARELLI JA, BORST LB, et al. Immune dysregulation and osteosarcoma: Staphylococcus aureus downregulates TGF-β and heightens the inflammatory signature in human and canine macrophages suppressed by osteosarcoma[J]. Vet Comp Oncol, 2020, 18(1): 64-75. DOI:10.1111/vco.12529 |
[14] |
XU HR, CHEN JJ, SHEN JM, et al. TYRO protein tyrosine kinase-binding protein predicts favorable overall survival in osteosarcoma and correlates with antitumor immunity[J]. Medicine, 2022, 101(39): e30878. DOI:10.1097/MD.0000000000030878 |
[15] |
DAI J, XU LJ, HAN GD, et al. Down-regulation of long non-coding RNA ITGB2-AS1 inhibits osteosarcoma proliferation and metastasis by repressing Wnt/β-catenin signalling and predicts favourable prognosis[J]. Artif Cells Nanomed Biotechnol, 2018, 46(sup3): 783-790. DOI:10.1080/21691401.2018.1511576 |
[16] |
CHEN LH, LIU JF, LU Y, et al. Complement C1q (C1qA, C1qB, and C1qC) may be a potential prognostic factor and an index of tumor microenvironment remodeling in osteosarcoma[J]. Front Oncol, 2021, 11: 642144. DOI:10.3389/fonc.2021.642144 |
[17] |
LI WH, DING ZY, WANG D, et al. Ten-gene signature reveals the significance of clinical prognosis and immuno-correlation of osteosarcoma and study on novel skeleton inhibitors regarding MMP9[J]. Cancer Cell Int, 2021, 21(1): 377. DOI:10.1186/s12935-021-02041-4 |