基于生物信息学预测甲状腺癌预后相关高风险糖酵解基因及其与患者预后的关系▲

张志勇 唐爱华 李双蕾 许淑华 高 美

(1 广西中医药大学研究生学院，南宁市 530000，电子邮箱：1916633986@qq.com；2 广西中医药大学第一附属医院内分泌科，南宁市 530000)

甲状腺癌是一种较常见的恶性肿瘤，近年来该病在全球范围内的发病率急剧上升[1]。甲状腺癌主要有甲状腺乳头状癌、甲状腺滤泡癌、甲状腺髓样癌、甲状腺间变性癌4种类型，其中甲状腺乳头状癌最为常见，占甲状腺癌的90%以上[2]。甲状腺癌的预后与其组织学类型密切相关，如甲状腺乳头状癌的10年总生存率约为98%，而甲状腺间变性癌的中位生存时间只有3～5个月[3]。对于甲状腺癌，手术切除是最重要的治疗方法，而不同的辅助治疗对某些病理类型的甲状腺癌是有效的，例如放射性碘治疗分化型甲状腺癌，以及化疗治疗甲状腺间变性癌。虽然大多数甲状腺癌患者预后良好，但约10%分化良好的甲状腺癌患者对放射性碘治疗无反应，而肿瘤分化程度差或未分化的患者更容易出现复发和死亡[4-5]。因此，寻找治疗甲状腺癌的新方法非常必要。糖酵解是肿瘤细胞获取能量的重要途径之一，阻断糖酵解途径有望成为治疗肿瘤的潜在方法。故本研究通过生物信息学技术预测甲状腺癌预后相关糖酵解高风险基因及其与患者预后的关系。

1 材料与方法

1.1 转录组数据及临床数据的下载 从肿瘤基因组图谱(The Cancer Genome Atlas,TCGA)数据库(https://portal.gdc.cancer.gov)下载甲状腺癌相关的转录组数据，即甲状腺癌相关基因及基因在临床样本(包括癌旁组织样本和肿瘤样本)中的表达量；对所获数据进行ID转换，即将样本名转换为基因名，获取每个样本中所包含的基因。同时下载患者的临床数据，包括患者的生存时间、生存状态(生存或死亡)、年龄、性别、临床分期等基本信息，剔除样本中信息不全的患者，得到最终的临床样本，获取其基因表达数据集文件及基因表型数据文件。

1.2 基因集富集分析 在基因集富集分析(Gene Set Enrichment Analysis，GSEA)数据库(http://www.gsea-msigdb.org/gsea/index.jsp)中搜集所有与糖酵解相关的基因集，将获取的基因表达数据集文件、基因表型数据文件及基因集文件导入GSEA软件，进行GSEA，筛选出P<0.05的基因集[6]。

1.3 甲状腺癌预后相关糖酵解基因的获取 应用Perl程序语言提取癌旁组织样本和肿瘤样本中糖酵解基因集(即1.2中筛选出的P<0.05的基因集)的表达量，采用Wilcoxon检验验证癌旁组织样本和肿瘤样本上述基因的表达量是否具有差异，筛选出P<0.05的差异性表达基因，输出基因表达量文件。将差异性表达基因的表达量与生存数据(生存状态、生存时间)合并，采用R软件对基因表达数据和生存数据进行相关性分析，得到与甲状腺癌预后相关的糖酵解基因。

1.4 网络模型的构建 根据甲状腺癌预后相关的糖酵解基因，采用R软件构建风险函数相关的Cox回归模型，得到糖酵解预后模型，由R软件自动计算出每个临床样本的风险值及所有临床样本风险值的中位值，当样本风险值高于中位值即判断为高风险样本，以此将患者分为高、低风险组，采用R软件绘制出两组患者的生存曲线图。并通过Cox回归模型分析结果获得各基因的风险比，当风险比>1时认为对应的基因即为高风险基因。

1.5 受试者工作特征曲线图、热图、时间-生存散点图的绘制及预后分析 采用R软件绘制出糖酵解预后模型预测甲状腺癌患者预后的受试者工作特征 (receiver operating characteristic，ROC)曲线图，当曲线下面积>0.7时，表示糖酵解预后模型预测的准确性高。根据糖酵解预后模型预测得到与甲状腺癌预后相关的基因，绘制热图及所有样本的时间-生存散点图，并进行预后分析。

2 结 果

2.1 临床样本的基本信息 共获得504例甲状腺癌临床样本，剔除3例信息不全的临床样本，共获得501例临床样本及56 753个基因数据。患者的生存时间为0～5 423 d，平均生存时间为1 031 d，中位生存时间为723 d；存活患者487例，死亡患者14例；发病年龄为15～89岁，平均发病年龄为45岁，中位发病年龄为46岁；女性患者366例，男性患者135例；临床分期Ⅰ期285例，Ⅱ期52例，Ⅲ期111例，Ⅳ期53例。

2.2 GSEA结果 通过GSEA数据库共搜索到5个与糖酵解相关的基因集，基因集名称分别为BIOCARTA_GLYCOLYSIS_PATHWAY、HALLMARK_GLYCOLYSIS、KEGG_GLYCOLYSIS_GLUCONEOGENESIS、REACTOME_GLYCOLYSIS、REACTOME_REGULATION_OF_GLYCOLYSIS_BY_FRUCTOSE_2_6_BISPHOSPHATE_METABOLISM。对这些基因集进行GSEA，结果显示只有HALLMARK_GLYCOLYSIS基因集的P值和校正P值小于0.05(见表1，仅列出该基因集的GSEA结果)，且富集分析图显示该基因在肿瘤样本中表达活跃(见图1)。

表1 GSEA结果

图1 HALLMARK_GLYCOLYSIS基因集的富集分析图

2.3 与甲状腺癌预后相关的糖酵解基因 获得147个差异性表达基因，将其表达量文件与生存数据合并后进行相关分析，共获得与甲状腺癌预后相关的糖酵解基因48个。

2.4 网络模型构建结果 将与甲状腺癌预后相关的48个糖酵解基因纳入Cox回归模型，构建糖酵解预后模型，得到与甲状腺癌预后相关的基因共15个，其中高风险基因共8个，分别为转化生长因子β诱导基因(transforming growth factor beta-induced gene,TGFBI)、溶质载体家族16成员3(solute carrier family 16 member 3,SLC16A3)、血小板型磷酸果糖激酶(phosphofructokinase platelet, PFKP)、丝氨酸/苏氨酸蛋白磷酸酶2A催化亚基β(serine/threonine-protein phosphatase 2A catalytic subunit beta，PPP2CB)、丙酮酸激酶M(pyruvate kinase M, PKM)、DNA损伤诱导转录子4(DNA damage inducible transcript 4,DDIT4)、碳水化合物磺基转移酶6(carbohydrate sulfotransferase 6，CHST6)、Quiescin巯基氧化酶1(Quiescin sulfhydryl oxidase 1，QSOX1)，见表2。根据临床样本的风险值及其中位值，将患者划分为高风险组、低风险组，两组患者的生存曲线见图2。由图可知：在疾病的早期，随着时间的延长患者的生存率降低，但在疾病后期，随着时间的延长患者的生存率逐渐趋于稳定，且高风险组的总体生存率低于低风险组(P=0.007)。

表2 甲状腺癌预后相关的高风险基因

图2 生存曲线图

2.5 ROC曲线分析结果 ROC曲线分析结果显示糖酵解预后模型预测甲状腺癌患者预后的曲线下面积为0.783，提示根据糖酵解预后模型预测甲状腺癌患者预后的准确率较高。见图3。

图3 ROC曲线

2.6 热图及时间-生存散点图 根据糖酵解预后模型得到的15个基因绘制出热图(见图4)和所有样本的时间-生存散点图(见图5)。从图4可以看出，表达量最高的4个基因分别是PKM、PPP2CB、PFKP、DDIT4，从图5可以看出，随着根据Cox回归模型计算出的风险值的增加，死亡的患者数增多。

图4 热图

图5 时间-生存散点图

3 讨 论

既往研究显示，肿瘤细胞比正常细胞需要更多的葡萄糖，即使在氧气充足的条件下，肿瘤细胞也是通过糖酵解进行葡萄糖代谢，而不是通过线粒体氧化磷酸化产生腺嘌呤核苷三磷酸，这就是“瓦伯格效应”[7]。然而，肿瘤细胞特殊的代谢模式不是被动的改变，而是遗传表达的积极改变，其改变了肿瘤发生和侵袭性的代谢模式。有研究表明，甲状腺癌的发生与葡萄糖摄取增加有关[8]。肿瘤细胞膜上的单羧酸转运蛋白(monocarboxylate transporter，MCT)与甲状腺癌的预后有关[9]，提示肿瘤代谢的关键分子可能是影响甲状腺癌预后的因素。

本研究通过生物信息学分析预测出与甲状腺癌预后相关的糖酵解高风险基因共8个，从ROC曲线可以看出，通过Cox回归构建的糖酵解预后模型预测甲状腺癌患者预后的准确性较高。研究显示，缺氧条件下源自肿瘤的外泌体包含信息性miRNA，其参与癌细胞与癌旁细胞的相互作用，从而促进肿瘤微环境的组织重塑[10]。与正常甲状腺滤泡细胞相比，从缺氧的甲状腺乳头状癌、BCPAP细胞及KTC-1细胞中分离出的外泌体增强了人脐静脉内皮细胞的血管生成；其次，在缺氧条件下，甲状腺乳头状癌BCPAP细胞的外泌体中miR-21-5p表达上调，而miR-21-5p可直接靶向并抑制TGFBI基因的表达，从而增加内皮细胞血管的形成，从敲除miR-21-5p的低氧BCPAP细胞中分离出的外泌体促进血管生成的作用减弱，提示缺氧的甲状腺乳头状癌细胞可通过外泌体miR-21-5p/TGFBI途径来促进血管生成[11]。国外学者发现，与滤泡腺瘤相比，TGFBI在滤泡癌中过表达[12]。PKM基因广泛分布于多种组织中，有M1和M2两种亚型。在肿瘤细胞中PKM2被PKM1替代以逆转warburg效应，从而减少乳酸的产生和增加氧气的消耗[13]。有学者发现，人甲状腺癌BCPAP细胞和TPC1细胞中PKM2 mRNA的表达水平高于非肿瘤细胞[14]。Bikas等[15]也发现PKM2在甲状腺癌FTC-133细胞和BCPAP细胞中过表达，且呈糖酵解依赖性。毒理学研究表明，DDIT4基因在甲状腺癌中非常活跃[16]。SLC16基因家族有14个成员，其中SLC16A1、SLC16A3、SLC16A7和SLC16A8这4个基因分别编码MCT1、MCT4、MCT2和MCT3，参与各种组织中L-乳酸、丙酮酸、短链脂肪酸和单羧酸盐药物的转运[17]，使用MCT1阻断剂能特异性靶向甲状腺癌ATC细胞，阻断癌细胞对丙酮酸的摄取[18-19]。但关于SLC16A3、PFKP、PPP2CB、CHPF2、QSOX1基因与甲状腺癌的相关性，目前尚未见文献报告，需进一步研究探讨。

肿瘤的糖酵解代谢能力对患者的预后有重要意义。糖酵解旺盛的肿瘤生长快、侵袭性强、复发迅速，患者的预后往往更差[20]。本研究的生存曲线显示，高风险组的生存率较低风险组下降，表明携带高风险基因的甲状腺癌患者的预后较差。国外有研究显示，MCT1可作为甲状腺癌预后的独立预测指标，抑制MCT1表达可使肿瘤中的恶性低氧细胞“间接饥饿”[19]。然而影响肿瘤预后的因素错综复杂，糖酵解代谢只是其中的一个环节，故通过糖酵解基因来判断患者的预后并不一定完全正确，抑制糖酵解高风险基因的表达，可能有助于改善甲状腺癌患者的预后，但仍需大量实验研究验证。

综上所述，甲状腺癌糖酵解途径与TGFBI、PKM、DDIT4、SLC16等高风险基因密切相关，抑制高风险基因的表达，可能有助于改善甲状腺癌患者的预后。但本研究仅在数据库的基础上挖掘甲状腺癌糖酵解途径的基因，有待今后进一步临床实验研究验证结果的可靠性。