甲状腺癌是一种常见的肿瘤，可发生于所有年龄段^[1]。目前，针对甲状腺癌的诊断主要依靠超声和经皮细针穿刺，治疗手段上则以手术和I¹³¹放射治疗为主。虽然，这种成熟的诊断和治疗方案得到了广泛的认可，但是，关于甲状腺癌的预后分析仍然是学者们探讨的热点。

众所周知，肿瘤细胞内微环境的紊乱多由供氧不足导致。缺氧可使基因和蛋白质的表达发生异常，也可能在遗传不稳定性、肿瘤发生和进展中发挥作用。已有研究^[2]表明，缺氧引起肿瘤微环境的变化，促进炎症发生、免疫抑制和治疗抵抗，从而诱发肺癌。

近年来，研究者们致力于利用缺氧相关标志物评估肿瘤的预后。例如，3个缺氧相关基因（PDSS1、CDCA8和SLC7A11）被用于构建肝细胞癌的预后、复发和诊断模型^[3]。基于缺氧相关基因，开发了缺氧风险模型来评估胶质瘤患者的预后^[4]。这些缺氧相关预后独立模型的建立对癌症患者的综合治疗做出了重大贡献。本研究旨在为甲状腺癌患者的综合治疗提供新的、可靠的预后标志物。

1 材料与方法 1.1 缺氧相关基因的获得

甲状腺癌样本的转录组数据和临床数据从肿瘤与癌症基因组图谱（the cancer genome atlas，TCGA）网站下载，其中，包括510份甲状腺癌样本数据和58份癌旁样本数据，以供进一步分析。根据错误发现率（false discovery rate，FDR） < 0.05和| log2FC | > 1，获得2 215个差异表达基因，结合从hypoxiaDB数据库（此数据库为人类缺氧相关调节蛋白的研究提供了一个完整的最新数据库）下载的差异表达基因，共获得了373个差异表达的缺氧相关基因。

1.2 功能分析

为进一步了解基因本身的功能，以及基因可能参与的通路，采用clusterProfiler进行基因本体（gene ontology，GO）分析及京都基因和基因组百科全书（Kyoto encyclopedia of genes and genomes，KEGG）分析，以研究差异表达的缺氧相关基因之间的功能相关性。

1.3 蛋白质-蛋白质相互作用网络的构建

将缺氧相关基因映射至STRING数据库，以创建一个交互网络来阐明差异表达基因之间的关联。随后，采用Cytoscape软件对蛋白质-蛋白质相互作用网络进行了可视化处理。

1.4 缺氧相关标志物预测模型的构建与验证

用Cox单因素回归分析鉴定与总生存时间相关的差异表达缺氧相关基因，并结合每个基因特定回归系数（β），构建患者的风险评分公式。风险评分= β（基因1）×表达量（基因1）+ β（基因2）×表达量（基因2）+……+β（基因n）×表达量（基因n）。

从整个甲状腺癌TCGA数据库中随机选择50%甲状腺癌样本（n = 255）作为实验组，其余50%甲状腺癌样本（n = 255）作为验证组，在2组根据风险评分公式，以中位风险评分为临界点，将患者分为低危组和高危组。用Kaplan-Meier分析评估2组的生存差异，并用log-rank统计方法进行比较。受试者操作特征（receiver operating characteristic，ROC）曲线和Harrell一致性指数用于评估模型预测的准确性。最后，采用单变量Cox回归分析来确定风险评分、年龄、肿瘤TNM分期以及性别之间的独立预后功能。

1.5 预测模型的内部验证

为了进一步验证该模型的预测能力，从整个甲状腺癌TCGA数据库中随机选择了50%甲状腺癌样本作为内部验证数据集（n = 255）。采用C指数评估所建立模型的性能。

1.6 统计学分析

采用R语言分析TCGA数据库。绘制Kaplan-Meier生存曲线，并通过log-rank检验进行比较。多因素Cox回归分析用于建立缺氧相关基因的预后预测模型。所有统计分析均使用R语言进行。所有统计检验均为双侧检验，P < 0.05为差异有统计学意义。

2 结果 2.1 差异表达缺氧相关基因的鉴定

从TCGA数据库中获得2 215个差异表达的甲状腺癌基因，结合hypoxiaDB数据库，共获得373个差异表达的缺氧相关基因。图 1A、1B所示为373个差异表达的缺氧相关基因中的前30个基因。其中，只有1个上调基因（CCND1），而下调基因有29个。采用箱氏图来总结甲状腺癌和癌旁组织之间30个差异表达的缺氧相关基因的表达模式，见图 1C。

2.2 差异表达的缺氧相关基因的功能富集分析

对差异表达的缺氧相关基因进行功能富集分析有助于理解其生物学作用。GO富集分析的前10个生物学过程和KEGG富集分析的前10个通路分析如图 2所示。差异基因的生物学过程主要涉及“细胞外结构组织”“含有胶原的细胞外基质”和“细胞外基质结构成分”（图 2A）。而差异基因的通路分析主要涉及“PI3K-Akt信号通路”“细胞因子-细胞因子-受体相互作用”和“MAPK信号通路”（图 2B）。

2.3 蛋白质-蛋白质相互作用网络的构建与重要基因模块

为了充分了解缺氧相关基因的差异表达，利用Cytoscape软件构建了一个交互式蛋白质-蛋白质相互作用网络。采用cytoscape软件中的分子复合物检测（MCODE）工具来识别重要的基因模块（hub genes），并筛选出3个重要的基因模块。

2.4 缺氧相关基因的鉴定与预后

采用单变量Cox回归分析确定与甲状腺癌患者预后相关的缺氧相关基因。森林图显示有23个基因在甲状腺癌患者中具有显著的预后价值（P < 0.05，图 3A）。进一步的多变量Cox回归分析共鉴定出11个基因（P < 0.05，图 3B）。

2.5 多因素Cox回归分析预测公式的建立

根据以下公式建立基因风险评分：遗传风险评分=0.792 6×SLC6A8+0.806 4×STC1+0.7867×PTGIS-1.057×SLC24A3-0.349 2×DPP4-0.749 2×ANK2-1.204 0×APOE+0.521 1×PIM1+3.085 0×PHKG1+0.547 3×AKR1C3+1.171 4×NDRG1。根据缺氧相关基因的预后公式，在实验组和验证组中分别确定了这些基因在不同风险人群中的分布以及患者的生存率。为了确定11个缺氧相关基因在预测甲状腺癌患者临床预后中的作用，进一步绘制了Kaplan-Meier生存曲线，以分析实验组和对照组中高危组和低危组的不同生存时间。Kaplan-Meier分析结果显示，高危组患者的生存率明显低于低危组，见图 4。

采用11个与缺氧相关的基因构建了2组甲状腺癌患者1年、3年和5年生存率的ROC曲线，来评估模型的预测性能（图 5）。在实验组中，ROC曲线下面积（area under curve，AUC）分别为0.943（1年）、0.897（3年）和0.831（5年）。而在验证组中，AUC分别为0.944（1年）、0.964（3年）和0.992（5年）。

对甲状腺癌患者的临床数据（包括风险评分）的预后价值进行单因素Cox回归分析（表 1），结果显示，构建的风险评分在预后评估中具有较高的预测价值。

3 讨论

近年来，甲状腺癌的发病率呈上升趋势，对人类健康构成了重大威胁。在美国，男性和女性甲状腺癌的检出率均有所上升，从1975年的每10万人中可检出4.9例上升至2014年的每10万人中14.3例^[4]。同时，由于5年特异性生存率差的低分化甲状腺癌的存在^[5]，对甲状腺癌的预后研究更为关键。

随着对缺氧引起的肿瘤微环境变化的深入了解以及TCGA和hypoxiaDB数据库的完善，学者们将更多的目光放在了通过缺氧相关基因来预测肿瘤的预后^[6-7]，近来新兴的证据都在强调免疫及炎症微环境在甲状腺癌中的重要性^[8-9]。然而，还没有一项研究致力于通过缺氧相关基因来对甲状腺癌患者进行预后分析。

目前，对甲状腺癌缺氧微环境的研究大多集中在单个基因上^[10-12]。本研究中鉴定的11个缺氧相关基因中，部分基因已被证实与肿瘤的预后密切相关。研究^[13]显示，CrT（SLC6A8）基因的敲除导致肌酸摄取不足，同时损害抗肿瘤T细胞的免疫功能。补充肌酸可显著抑制肿瘤生长，并与PD-1/PD-L1阻滞剂协同抑制肿瘤生长^[14]。载脂蛋白E作为高脂血症的主要病因，被证明在高脂血症模型中，可通过运动改善肿瘤的缺氧微环境，同时减缓原发性和继发性EO771乳腺肿瘤的形成^[15]。

此外，PTGIS和DPP4在生物信息学研究中也得到了证实。在先前与肝细胞肝癌相关的新生物标志物研究中，PTGIS是通过mRNA表达网络分析确定的21个核心基因之一，可能是抑制肝细胞肝癌的潜在治疗靶点^[16]，而DPP4也被证明是NF-κB 15的靶基因^[17]。

STC1、PIM1和NDRG1这3种基因在癌症中的研究更加成熟。由于STC1可能通过Bcl-2的失调在缺氧性胃癌中发挥致癌作用，因此被视为可能的胃癌潜在治疗靶点^[18]。在缺氧微环境下，miR-124和miR-144的下调可通过减弱PIM1的抑制作用使细胞发生缺氧，而增加患前列腺癌的风险^[19]。NDRG1作为调节脂质代谢的关键基因可促进乳腺癌的侵袭性。因NDRG1与乳腺癌预后不良之间的关系密切，故NDRG1成为了乳腺癌中一种有前途的治疗靶点^[20]。

总之，本研究基于对差异表达的缺氧相关基因表达谱和相应临床特征的综合分析，确定了多种甲状腺癌的预后标志物。并通过构建一个新的风险评分模型，以有效地评估甲状腺癌患者的预后。本研究的局限性在于其为回顾性研究，今后应该进行更多的前瞻性研究，以验证缺氧相关基因的预后功能。当然，这也需要多中心数据来证实本研究的发现。