基于骨肉瘤核心驱动基因筛选的生物信息学分析和患者生存期预测基因模型的构建

李苇航, 丁子毅, 王 栋, 潘益凯, 刘玉辉, 张世磊, 李 靖, 闫 铭

（1.中国人民解放军空军军医大学西京医院骨科，陕西 西安 710032；2.中国人民解放军空军军医大学航空航天医学系航空航天医学训练教研室，陕西 西安 710032；3.中国人民解放军空军军医大学航空航天医学系航空航天临床医学中心 教育部航空航天医学重点实验室，陕西 西安 710032）

目前，对于骨肉瘤（osteosarcoma，OS）的临床治疗主要包括手术、放疗和化疗。随着新辅助化疗和保肢手术的引入及发明，OS 患者的总体生存率有了显著提高。尽管手术和靶向化疗取得了进展，但由于OS 高度侵袭性和进展迅速的特点，约20% 的患者在诊断时已发生远处转移［1-2］。肿瘤术后复发或转移仍然是OS 患者长期生存率低的主要原因，其总体预后不良。研究［3］显示：OS 的发病机制与基因的表达密切相关。而关于OS 致病基因集的探索及其与患者临床预后的关联性研究尚欠缺。在临床的早期阶段，从分子水平探讨OS 发生发展的靶基因以及治疗靶点有助于临床工作者对该病进行早发现、早诊断和早治疗。进一步探讨用于OS 早期诊断和预后评估的生物标记物是研究者亟需解决的问题。

近年来，随着高通量测序技术和基因微阵列芯片技术的快速发展，生物信息学在识别肿瘤的驱动基因、阐明致癌机制及预测预后等方面展现出了广阔的应用前景［4-5］。本研究采用生物信息学分析方法，首先从基因表达汇编（Gene Expression Omnibus，GEO）数据库查找不同样本中OS 基因的表达谱芯片数据，通过对比分析的统计方法，获得OS 的差异表达基因（differentially expressed genes，DEGs），随后对DEGs进行基因本体论（Gene Otology，GO）和京都基因和基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）信号通路分析，从分子功能、生物学过程和细胞成分方面揭示其发病机制。构建蛋白-蛋白相互作用（protein-protein interaction，PPI）网络，探讨导致癌变的相关蛋白之间的关系。此外，本研究在肿瘤基因组图谱（The Cancer Genome Atlas，TCGA）数据库中获得OS 患者的临床信息和转录组数据，通过Kaplan-Meier（K-M）生存分析全面揭示驱动基因、不同种族和不同性别OS 患者的预后情况，探讨OS 发生发展的分子机制及其潜在的治疗靶点以及可能影响预后的相关因素，为临床工作者提供参考。

1 资料与方法

1.1 基因芯片获取和数据预处理于美国国立生物技术信息中心（NCBI）的GEO数据库（https：//www.ncbi.nlm.nih.gov/geo/）中获取OS芯片数据GSE12865、GSE14359 和GSE36001，分别基于GPL6244、GPL96 和GPL6102 平台获得探针注释信息。原始数据通过R语言“affy”和“affyPLM”包分析［6］，将原始的CEL 文 件通过RMA 算法（“rma” 函数）对数据进行背景矫正、质量控制及标准化处理，转换成探针表达矩阵。随后将3个GSE数据集中的探针基于各自的注释平台转换成相应的基因名称，移除无对应基因名的探针集，得到3个数据集的基因表达矩阵。

1.2 临床患者的数据收集通过R 语言的“TCGAbiolinks”包［7］下在TCGA数据库（https：//portal.gdc.cancer.gov/）中获得OS 患者的所有相关临床信息，包括患者年龄、性别、种族、末次随访日期和生存状态等因素。共收集了379个样本相关的临床记录信息，其中88个样本有相应的基因表达谱。

1.3 OS 的DEGs筛选采用R 语言的“limma”包对骨肉瘤样本和正常样本组织的基因表达谱进行分析，同时满足|logFC|（fold change）＞1 且P＜0.05 的被鉴定为DEGs，认为差异有统计学意义。将3个数据集分别获得的DEGs 进行韦恩图分析（R语言“VennDiagram”包），得到更加精确的共同DEGs的结果。随后通过“gplots”包中的“heat.map2”函数绘制3个数据集的热图，直观地显示出共同的DEGs在不同组织样本中的表达情况。

1.4 OS驱动基因的分子功能和通路富集分析采用DAVID（http：//david.abcc.ncifcrf.gov/）基因功能注释在线数据库工具［8］，对筛选出的DEGs进行主要功能和信号通路分析，探讨其DEGs 背后潜在的真正生物学意义。GO分析研究DEGs 富集的具体分子功能、生物学过程和细胞成分。KEGG是基因组信号通路关联的基础，用于分析有关信号通路改变的信息。在DAVID 数据库中将DEGs 进行上述分析，以P＜0.05 作为差异有统计学意义的阈值。基因集富集分析（Gene Set Enrichment Analysis，GSEA）进一步探讨DEGs 主要作用的细胞信号通路。

1.5 PPI 网络分析STRING 数据库（https：//string-db.org/）是一个预测不同基因表达的蛋白与蛋白之间相互作用的在线检索工具［9］，既包括不同蛋白之间直接的物理相互作用，也包括了蛋白之间的功能相关性。采用STRING 数据库对筛选出的DEGs 进行PPI分析，构建相互作用网络，随后将构建好的网络导入Cytoscape 3.8.0 软件中，采用“cytohubba” 插件中的MCC 算法进一步筛选出关系最紧密的靶基因。“cytohubba” 插件中的MCC算法旨在从复杂庞大的相互作用网络中发现最相关的关键靶点及子网络，这是识别鉴定靶基因的最有效的方法之一［10］。

1.6 OS 患者临床样本的生存分析采用R 语言下载TCGA 数据库中的OS 临床样本，共88个样本有相应的基因表达谱。当前生存状态和随访时间资料完整的患者被纳入本研究，共纳入86例患者。患者基于关键基因表达值的中位数分为高表达和低表达两组，采用R语言的“survival”和“survminer” 包进行K-M 生存分析全面探讨不同驱动基因表达水平、不同种族和不同性别OS 患者的预后情况。

1.7 OS 患者总生存期的临床预测应用采用单因素和多因素Cox 回归分析评估驱动基因集的预后能力，基于多因素Cox 回归分析的结果绘制列线图。对驱动基因进行建模来预测OS 患者预后（3年生存率、5年生存率和10年生存率），评估其预测能力，并绘制校准图以验证其预测准确性。采用C 指数（Concordance index）描述建模的辨识度，C 指数值越高则建模预测越可靠。

2 结 果

2.1 数据处理和DEGs 的鉴别采用R 语言对CEL 文件进行背景矫正和标准化处理，标准化后的基因表达谱以小提琴图呈现。见图1。

图1 经过标准化校正后样本表达量的小提琴图Fig.1 Violin diagrams of expressions of samples after standardized correlation

基于|logFC|（fold change）＞1 且P＜0.05 的截断值，GSE12865 共鉴别出1 632个DEGs，其中有605个上调基因和1 027个下调基因；GSE14359鉴别出766个DEGs，其中有245个上调基因和521个下调基因；GSE36001共鉴别出358个DEGs，其中有203个上调基因和155个下调基因。通过韦恩图分析，对这3个数据集的DEGs 取交集以获得更加精确可靠的结果（图2A），共有18个基因鉴定为共同DEGs。层次聚类分析结果显示：2组不同样本（肿瘤组与正常组）之间的共同DEGs 的表达量差异均有统计学意义（图2B～2D）。共同上调基因有CD93、ARHGDIB、LYZ、CD74、LAPTM5、FCER1G、HLA-DPA1、HCLS1、GIMAP4、RNASE1、SNRPA1 和TYROBP；共同下调的基因有LTBP2、PNMA2、DDAH1、TAGLN、TLN1 和ANPEP。

图2 DEGs 的筛选与不同样本中DEGs 的表达水平Fig.2 Screening of DEGs and expression level of DEGs in different samples

2.2 DEGs 的功能和通路富集分析将3个微阵列芯片数据共同上调和下调的基因导入DAVID 在线工具中，进行了GO 和KEGG 富集分析。GO分析结果表明： 上调的DEGs 主要涉及细胞外基质分解、软骨内成骨、特定蛋白结合、成骨细胞分化、胶原分解代谢和发挥调节MHC Ⅱ类受体活性，参与构成细胞外基质及MHC Ⅱ类蛋白复合体，富集于细胞表面；KEGG 信号通路分析结果显示：上调的DEGs 主要参与多糖在癌症中的表达、细胞周期中DNA 的复制、PI3K-AKT 经典信号通路（图3A）；下调的GO和KEGG通路富集结果见图3B。GSEA分析结果表明：DEGs 在Notch 信号通路中起到明显作用，在OS 的发生发展中产生重要影响（图3C）。

图3 DEGs 的GO、KEGG 和GSEA分析结果Fig.3 Analysis results of GO，KEGG ，and GSEA of DEGs

2.3 PPI的构建和核心驱动基因的筛选本研究将筛选出的共同DEGs 输入STRING工具，并将所得到的基因作用网络导入Cytoscape 3.8.0 软件行进一步分析，利用 “Cytohubba”插件中的MCC 算法找出基因网络中的子网络，即相互作用最紧密的基因集，按照MCC方法计算得到的前10名基因分别 是TYROBP、LAPTM5、FCER1G、CD74、HCLS1、ARHGDIB、HLA-DPA1、CD93、GIMAP4 和LYZ。见图4，其中颜色越深表明该基因在OS 的发生发展中所起的作用越大。

图4 OS 的DEGs 之间相互作用关系Fig.4 Interaction relationships between DEGs of OS

2.4 OS 患者的生存曲线分析为了进一步验证本研究筛选出的驱动基因对OS 患者预后的影响程度，本研究采用PPI分析出的10个核心驱动基因为研究对象，对86例OS 患者进行K-M 生存曲线分析，基于基因中位表达值分为高表达组和低表达组的结果显示： ARHGDIB、CD74、FCER1G、HCLS1、HLA-DPA1 和TYROBP 基因低表达组与高表达组患者生存率比较差异有统计学意义（P＜0.05），提示这些基因与OS 患者的预后具有高度相关性，上述基因低表达组OS 患者预后良好，生存率高，生存期更长（P＜0.05），OS 患者的生存曲线见图5。

图5 OS 患者DEGs 的K-M 生存分析曲线Fig.5 K-M survival analysis curves of DEGs of OS patients

不同性别和不同种族OS 患者的生存曲线见图6，性别和种族均不会对OS 患者的预后产生影响（P＞0.05）。

图6 不同性别和种族OS 患者的K-M 生存分析曲线Fig.6 K-M survival analysis curves of OS patients with different genders and races

2.5 OS 患者总生存期的预测为了预测OS 患者的总生存期，为OS 患者预后转归提供指导，本研究从生存分析中选择3年生存期差异有统计学意义的6个基因组成基因特征集进行建模，通过检测患者样本相关基因的表达量预测患者的3年、5年和10年生存率。在多因素Cox 回归分析绘制的列线图模型中，C指数（0.71）显示出该模型良好的鉴别度，列线图结果见图7A。此外，本研究还对其预测能力进行了分析，校准图显示出该模型预测OS患者3年、5 和10年生存率均有高度的准确性（图7B）。

图7 基于6个基因特征集构建的列线图模型和校准图Fig.7 Nomogram model and calibration plot constructed by 6-gene feature sets

3 讨 论

OS 是骨科最常见的恶性肿瘤之一，占所有类型骨癌的20%～40%，是儿童和青少年最常见的原发肉瘤，每年（0.3～0.4）/100 000 人患病［11-12］。OS 的特征是在未成熟的骨中直接形成类骨样组织、类骨质和梭形基质细胞。OS主要临床表现为肿胀和骨痛，主要发生于四肢长骨，特别是膝关节［13］，而体质量减轻、面色苍白、盗汗、发烧和食欲不振等全身性症状相对少见［14］。近年来，尽管OS手术术式和辅助化疗技术取得了较大进步，但OS患者的总体生存时间仍然较短，5年生存率约为20%，且肿瘤复发率高［15］。研究［16］显示：OS患者生存率改善不理想的主要原因可能与OS的转移有关。因此需要在临床早期阶段进行诊断，避免延误诊断而导致转移的发生。精确的基因诊断是进行靶向治疗和评估预后的关键，为了寻找潜在的药物治疗靶点及新的OS标志物，需要进一步探讨OS 的发病机制。目前关于OS 的研究主要集中在OS 的转移机制上，而对于造成OS 发生发展的分子靶点的相关研究还不够深入。因此更多地了解OS 的发病机制、探索OS 的治疗靶点和临床早期诊断OS 是亟需解决的问题。

近年来，基因芯片测序和二代测序技术的快速发展，以及生命科学领域新兴交叉学科生物信息学的崛起，测序技术结合生物信息学的方法被广泛应用于基因组学及转录组学上研究疾病的发生发展及转移方面［4，17-18］，可用于分析各种疾病背后的致病基因及其真正的生物学意义。本研究针对OS 患者的临床大样本数据集，采用标准化的统计方法，从GEO 数据库中对OS 患者基因芯片数据集的mRNA 表达情况进行分析，将3个芯片分析出的DEGs 进行韦恩图绘制，得到更加精确的共同DEGs，并进行GO、KEGG 和GSEA 通路富集分析，随后对共同DEGs 构建PPI 网络，加之以基因权重分析。将筛选出的权重基因进行生存分析进一步分析患者预后，同时对3年生存率差异有统计学意义患者的基因特征集建模，绘制列线图及校准图以验证基因特征集的预测预后能力。

经过韦恩图分析后的共同上调DEGs 有CD93、ARHGDIB、LYZ、CD74、LAPTM5、FCER1G、HLA-DPA1、HCLS1、GIMAP4、RNASE1、SNRPA1 和TYROBP。GO分析包括了生物过程（biological process，BP）、分子功能（molecular function MF）和细胞成分（cellular component，CC）3个部分，研究结果显示：DEGs 主要涉及有丝分裂的核分裂、细胞外基质分解、软骨内成骨、成骨细胞分化、胶原分解代谢和发挥调节MHC Ⅱ类受体活性，参与构成细胞外基质和MHC Ⅱ类蛋白复合体，富集于细胞表面；KEGG 信号通路分析结果显示：上调的DEGs 主要参与多糖在癌症中的表达、细胞周期中DNA 的复制和PI3K-AKT 经典信号通路。GSEA分析结果表明：在Notch 信号通路中，DEGs 起明显作用，以上结果均预示DEGs和这些信号通路在OS 的发生发展中起重要作用。

肿瘤的发生不是单个基因变化的结果，而是由许多基因和信号通路共同激活所造成的。本研究在STRING 数据库中构建PPI，根据MCC 算法选取排名前10位的Hub基因，分别是TYROBP、LAPTM5、FCER1G、CD74、HCLS1、ARHGDIB、LA-DPA1、CD93、GIMAP4 和LYZ，随后将排名前10名的基因进行生存分析以进一步验证肿瘤患者预后情况，最终分析出对患者预后具有意义的基因特征集（ARHGDIB、CD74、FCER1G、HCLS1、HLA-DPA1 和TYROBP），共6个基因，该6个基因被认为是导致OS 形成最重要的基因。其中CD74、HCLS1和HLA-DPA1在既往关于OS 的研究［19-21］中有过报道，且与本研究的结论一致，证实了这些基因在OS 的发生发展中起主导作用。ARHGDIB属于Rho鸟苷二磷酸解离抑制剂（RHO guanosine diphosphate dissociation inhibitors，ARHGDIs）家族，以前被认为是一种次要的组织相容性抗原，是一种参与许多细胞活动的细胞内GTP结合蛋白［22］。研究［23］显示：ARHGDIB只在造血组织中表达，主要是在淋巴细胞中。随着研究的进展，ARHGDIB被发现也在其他器官和组织中广泛表达。研究［24］显示：ARHGDIB 是一种侵袭性的人类癌症标志物。研究［25］显示： ARHGDIB具有致癌作用，例 如ARHGDIB表达的缺失抑制了胰腺癌的侵袭和迁移能 力。研究［22］显示：ARHGDIB 通过激活肝细胞癌中的AKT 来促进细胞的侵袭和增殖。本研究发现ARHGDIB 在导致OS 的进展中同样发挥了重要作用。FCER1G 位于染色体1q23.3上，其编码的蛋白质是高亲和力IgE 受体和白细胞介素3（interleukin-3，IL-3）受体复合物的组成成分，FCER1G 主要参与介导肥大细胞的过敏性炎症信号，选择性地介导嗜碱性粒细胞产生白细胞介素4（interleukin-4，IL4），启动T 淋巴细胞向效应T 淋巴辅助细胞2亚群的转移［26-27］。FCER1G是一种先天免疫基因，参与多种疾病的发生进展，例如湿疹、脑膜瘤、儿童白血病、肾透明细胞癌和急性髓系白血病等［28-31］，可能通过影响免疫相关途径改善预后。此外，本研究也证实了FCER1G 与OS 的进展有紧密关联。TYROBP，也称DAP12，位于染色 体19q13.1，是一种Ⅰ型跨膜蛋白［32］。TYROBP 表达在免疫细胞上，包括巨噬细胞、单核细胞和骨细胞［33-34］，被认为是激活信号转导的关键成分，这与本研究的结论一致，即TYROBP 在OS组织中高表达。一项关于阿尔兹海默病的研究［35］显示：降低TYROBP 的表达量对神经具有保护作用，代表着一个重要治疗或是预防的机会，结合本研究分析出的TYROBP 在OS 中高表达的结论，上述研究均证实了TYROBP 的低表达对机体有积极影响。本研究通过生存分析证实了OS 患者的性别和种族的差别对患者的生存时间无影响。

根据这6个基因特征集的表达量，本文作者建立了预测生存期的模型，C 指数以及校准图均说明了该模型具有良好的预测能力。作为临床诊断性生物标志物、治疗靶点和预后指标，在临床上通过检测该基因特征集各自的表达量即可以对OS 进行鉴别诊断及预测患者的生存时间，有利于OS 的早发现、早诊断和早治疗。

本研究通过生物信息学的分析方法，从不同基因芯片数据中，探讨导致OS 进展的基因之间的联系和影响OS 发生的遗传机制，为OS 的预防、诊断和治疗提供了理论基础。ARHGDIB、CD74、FCER1G、HCLS1、HLA-DPA1 和TYROBP 基因特征集在OS 发生发展中发挥主要作用，其中ARHGDIB、FCER1G 和TYROBP 首次被证实与OS 进展有关，为OS 的早期诊断、靶向治疗及预后情况提供了新思路。本研究虽然基于临床大样本数据进行分析讨论，但仍缺乏具体的细胞和动物实验进行验证，尚存在不足之处，未来需要进一步的体内外实验进行完善。