甲状腺癌铁死亡预后风险模型的构建及其潜在机制的生物信息学分析

杨仁义, 彭书旺, 王永恒, 董宇轩, 段姗姗

（1.湖南中医药大学第一附属医院胃肠甲状腺血管外科，湖南 长沙 410007；2.湖南中医药大学研究生院，湖南 长沙 410208）

国家癌症中心发布的《2022 年全国癌症报告》［1］显示：甲状腺癌（thyroid cancer，TC）新发病例为20.3 万人，其中女性新发病例（15.3 万人）明显高于男性（5.0 万人），是发病率最高的内分泌恶性肿瘤，也是发病率上升最快且位列第7 位的常见恶性肿瘤。95%的TC 为分化型TC，外科手术是目前首选的治疗方式，同时放射性碘治疗（radioactive iodine，RAI）和内分泌辅助治疗也是提高患者生存预后和生活质量的有效手段［2-3］，但治疗后远期复发率仍较高，针对免疫检查点或分子靶向治疗药物开发大多数仅停留在疗效观察阶段［4-5］，尚未取得分子机制层面的突破。因此，研究TC 潜在的分子机制，寻找预后生存的生物标志物和治疗靶点，构建有价值的预后评估模型是亟待解决的问题。铁死亡是一种主要依赖于铁介导的氧化损伤和随后的细胞膜损伤的调节性细胞死亡方式［6］，可通过外源性或转运蛋白依赖，以及内源性或酶调节2 条主要途径启动，同时受铁积累增加、自由基产生、脂肪酸供应和脂质过氧化共同调控［7］，在一定程度上抑制TC细胞的增殖、侵袭及转移。研 究［8-10］显 示：铁 死 亡 相 关 基 因（ferroptosisrelated genes，FRGs）是TC 预后分子及潜在治疗靶点，能精准预测TC 患者生存预后，有效指导TC 患者的个体化治疗。WANG 等［11］鉴定出醛酮还原酶家族1 成员 C3 （Aldo-Keto reductase family 1 member C3，AKR1C3）、BH3 相互作用域死亡激动剂（BH3 interacting domain death agonist，BID）、F-Box和 WD-40结构域蛋白7（F-Box and WD repeat domain containing 7，FBXW7）、谷胱甘肽过氧化物酶4 （glutathione peroxidase 4，GPX4）和丝裂原活化蛋白激酶激酶5（mitogen-activated protein kinase kinase kinase 5，MAP3K5）具有预后价值，同时敲低AKR1C3 可增强TC 细胞的增殖、侵袭和迁移能力，为TC 预后预测提供了新思路，但未对其相关分子机制进行深入探讨。本研究采用生物信息学方法，提取GeneCards 和FerrDb 数据库中FRGs，筛选TC 差异预后铁死亡相关基因（prognosis and ferroptosis related genes，PFRGs），构建风险评分及预后模型，评估TC患者生存预后，进行共表达及富集分析，研究TC的FRGs调控网络及作用机制。

1 资料与方法

1.1 TC 差 异 预 后 基 因（TC-differentiallyprognostic gene，TC-DPGs）获取数据来源：从TCGA 数据库（https：//portal.gdc.cancer.gov/）下载TCGA-THCA 队列中TC 转录组及临床数据。差异分析：采用R 软件DESeq2 包［12］，对TCGATHCA 转录组数据进行差异分析，以|log2FC|>1且校正P<0.05 为筛选条件，筛选TC 差异表达基因 （thyroid cancer differentially expressed genes，TC-DEGs）。预 后 分 析：采 用R 软件survminer 包和survival 包［13］，对TCGA-THCA 临 床 数 据 进 行预后分析，以TC 患者总生存（overall survival，OS）时间为预后参数，以Cox 回归P<0.05 为筛选条件进行单因素Cox 回归分析，筛选TC 预后基因（TC prognosis gene，TC-PGs）。通过韦恩图在线取交集工具（http：//bioinformatics.psb.ugent.be/webtools/Venn/）取“TC-DEGs”与“TCPGs”两者的交集，获取TC-DPGs。

1.2 TC 差异PFRGs 筛选数 据 来 源：通 过GeneCards 数 据 库 （https：//www.genecards.org/）［14］，以“Ferroptosis”为 检 索 词，同 时 从FerrDb 数 据 库 （http：//www.zhounan.org/ferrdb）［15］下载FRGs，整合2 个数据库基因，最终获取FRGs。通过韦恩图在线取交集工具，取“TC-DPGs”与“FRGs”交 集，获 取TC 差 异PFRGs。

1.3 TC 差异PFRGs 表达分析数据来源：从TCGA 数据库下载TCGA-THCA 队列中TC 转录组数据，从GTEx 数据库（https：//gtexportal.org/home/）下载正常甲状腺组织表达数据。将TCGA 数据库中癌旁组织样本与GTEx 数据库中正常甲状腺组织样本作为对照组，将TCGA 数据库中TC 组织样本作为肿瘤组，进行非配对样本秩和检验；对包含肿瘤组织与癌旁组织的患者样本进行配对样本秩和检验，比较肿瘤组与对照组TC 差异PFRGs ［CD44、膜 联 蛋 白A1 （Annexin A1，ANXA1）和 核 受 体 亚 家 族4A 类 成 员1 （nuclear receptor sub family 4 group A mumber 1，NR4A1）］表达的差异。为进一步分析TC 差异PFRGs（CD44、ANXA1 和NR4A1）蛋白在组织中表达情况，通过人类蛋白图谱（The Human Protein Atlas，HPA）数 据 库［16］（https：//www.proteinatlas.org/），检索CD44、ANXA1 和NR4A1蛋白的免疫组织化学切片图像，比较肿瘤组织与正常组织蛋白表达差异。

1.4 TC 差异PFRGs 生存预后分析采用R 软件timeROC 包，绘制时间依赖性受试者工作特征（time-receiver operating characteristic，time-ROC）曲线，分别评估TC 差异PFRGs CD44、ANXA1和NR4A1 对TC 患者的1、3 和5 年生存预测效能。当ROC 曲线下面积（area under curve，AUC）>0.5 时，AUC 越 接 近1，对TC 患 者 的1、3 和5 年生存预测效能越好，分子的表达促进死亡发生；当AUC<0.5时，AUC越接近0，对TC患者的1、3 和5年生存预测效能越好，分子的表达促进生存发生。采用R 软件survival 包和survminer 包，绘制Kaplan-Meier 曲线，根据基因表达的中位数分为高表达和低表达2 组，分别评估TC 差异PFRGs CD44、ANXA1 和NR4A1 对TC 患者生存状况的影响。

1.6 TC 差异PFRGs Nomogram 图构建及验证采用R 软件rms 包，使用多因素Cox 回归分析方法建立Nomogram 预后模型，以评估和预测TC患者OS 的概率［17］。综合分析“1.5”中单因素Cox 回归分析结果，纳入多因素Cox 回归分析中P<0.05 的临床特征作为因子，构建预测TC 患者1、3 和5 年OS 概率的Nomogram 预后模型。采用R 软件rms 包和survival 包，对Nomogram 预后模型进行Calibration 分析验证，以评估Nomogram 预后模型的预测能力。

1.7 TC 差异PFRGs 与蛋白共表达相关性分析

采 用R 软 件limma 包［18］，对TC 差 异PFRGs CD44、ANXA1 和NR4A1 进 行Spearman 相 关 性 分析，以相关系数R>0.7 且P<0.05 为筛选条件，筛选TC 差异PFRGs 的共表达基因，绘制共表达热图。整合CD44、ANXA1 和NR4A1 的共表达基因，通过STRING（https：//cn.string-db.org/）［19］数据库进行蛋白互作分析，采用Cytoscape 3.9.0 绘制蛋白-蛋白互作（protein-protein interaction，PPI）网络图，利用cytoHubba 插件［20］进行度中心性（degree）、接近中心性（closeness）和中介中心性（betweenness）拓扑分析，探寻TC 铁死亡核心网络。

1.8 TC 差 异 PFRGs 基 因 本 体 论（Gene Ontology，GO）和京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析采 用R 软 件clusterProfiler［21］、org.Hs.eg.db 包，对TC 差 异PFRGs CD44、ANXA1 和NR4A1 的共表达基因进行GO 和KEGG富集分析，以校正P<0.05 且q-value 的值<0.2 为筛选条件，聚焦TC 铁死亡相关的生物进程（biological process，BP）、细 胞 组 分 （cell component，CC）、分子功能（molecular function，MF）和KEGG 通路，并构建“KEGG-靶点”网络图。

1.9 统计学分析采用R 软件（3.6.3）中的Bioconductor 软件包进行数据处理和统计学分析。TCGA 数据库中转录组数据不符合正态分布，采用Wilcoxon检验评估肿瘤组和对照组之间PFRGs表达的差异，采用单因素和多因素Cox 回归分析评估PFRGs 表达和临床数据等与TC 患者OS 之间的相关性。以P<0.05 为差异有统计学意义。

2 结 果

2.1 TC 差异预后基因获取结果以|log2FC|>1 且校正P<0.05 为筛选条件，鉴 定 出TC-DEGs 共6 773 个，其中3 317个在TC中表达上调，3 456个在TC 中表达下调，绘制火山图（图1A）及差异排序图（图1B）。以P.cox<0.05 为筛选条件，鉴定出TC-PGs 共1 444 个，其 中HR>1 的 危 险 基 因755 个，HR<1 的保护基因有689 个。取“TCDEGs”与“TC-PGs”两者的交集，韦恩图可视化，获取TC-DPGs 共343 个（图1 C）。

图1 TC 差异预后基因筛选结果Fig.1 Screening results of TC differential prognosis genes

2.2 TC 差异PFRGs 筛选结果从GeneCards 数据库中获得FRGs 442 个，从FerrDb 数据库中获得FRGs 384 个，2 个 数 据 库 存 在167 个 相 同FRGs，保留唯一值后，共获取659 个FRGs。取“TCDPGs”与“FRGs”两者的交集，韦恩图可视化，获 取PFRGs 共3 个（图2），分 别 为CD44、ANXA1 和NR4A1，差异分析及预后分析结果见表1。TC 中CD44 和ANXA1 差异表达上调（|log2FC|>1）是保护性因素（HR<1），NR4A1 差异表达 下 调 （|log2FC| <－1）是 危 险 性 因 素（HR>1）。

表1 TC 差异PFRGs 差异和预后分析Tab.1 Differences and prognostic analysis on TC differential PFRGs

图2 TC 差异PFRGs 韦恩图Fig.2 Venn diagram of TC differential PFRGs

2.3 TC 差异PFRGs 表达情况配对或非配对样本t检验结果（图3A 和3B）显示：与正常组织比较，TC 差 异PFRGs CD44 和ANXA1 在TC 中 表达上调，NR4A1 在TC 中表达下调。同时HPA 数据库免疫组织化学图片结果（图4）显示：TC 组织样本中CD44 和ANXA1 蛋白过表达，NR4A1 蛋白低表达，此结果与TCGA 数据库结果相互验证。

图3 TC 差异PFRGs 的表达Fig.3 Expressions of TC differential PFRGs

图4 免疫组织化学法检测TC差异PFRGs蛋白表达（×400）Fig.4 Protein expressions of TC differential PFRGs genes detected by immunohistochemistry(×400)

2.4 TC 差异PFRGs 生存预后情况time-ROC 生存预测结果（图5A～5C）显示：TC 患者1、3 和5 年生存预测中，ANXA1 均表现出相对最佳的预测效能，其中5 年的预测效价（AUC=0.200、0.221 和0.198）最佳；CD44 预测5 年生存预后具有最好效能，基因表达在TC 患者早期（1 年）促进死亡发生，在中后期（3 和5 年）促进生存发生；ANXA1 表 达 在TC 患 者1、3 和5 年 中 均 促 进 生 存发生；NR4A1 表达在TC 患者1、3 和5 年中均促进死亡发生。

Kaplan-Meier 曲线生存状况结果（图5D～5F）显示：TC 差异PFRGs CD44、ANXA1 和NR4A1的表达与生存预后相关（P=0.048、0.005 和0.036），其中CD44 和ANXA1 的表达是保护性因素，NR4A1 的表达是危险因素，即低表达CD44 和ANXA1 与高表达NR4A1 具有较差的生存结局。

图5 TC 差 异PFRGs 的time-ROC 和Kaplan-Meier 曲 线Fig.5 Time-ROC and Kaplan-Meier curves of TC differential PFRGs

2.5 TC 差异PFRGs 独立预后分析结果对TC差异PFRGsCD44、ANXA1 和NR4A1 的表达进行多因素Cox 回归分析，计算出3 个基因比例风险回归模型的Coef 相关系数，即ANXA1（coef）=－0.361，NR4A1（coef）=0.316，CD44（coef）=0.008。将Coef 相关系数与3 个基因的表达值纳入风险评分计算公式，计算TC 患者的风险评分。根据患者风险评分中位数（－0.049），将TCGA-THCA 甲状腺患者分为高风险组255 例和低风险组255 例。

高和低风险组生存差异比较（图6），Log-rank检 验 结 果 显 示：HR=6.87，95%CI：2.58～18.30，P=0.003；Cox 回归分析结果显示：HR=6.88，95%CI：1.56～30.28，P=0.011，高风险组和低风险组患者生存预后比较差异有统计学意义，同时高风险组较低风险组具有更差的生存预后。time-ROC 显示出风险评分能良好地预测TC患者1、3 和5 年的生存预后，1、3 和5 年的AUC值分别为0.761、0.767 和0.722，其中预测3 年的生存预后效价最优。

图6 低风险组和高风险组TC 患者Kaplan-Meier 曲线（A)及time-ROC 曲线（B）Fig.6 Kaplan-Meier curve(A) and time-ROC curve(B) of TC patients in low risk group and high risk group

为进一步评估风险评分与其他临床特征的预后关系，采用单因素和多因素Cox 回归分析，以单因素Cox 回归分析的P<0.1 作为纳入多因素Cox 回归分析的条件，对临床特征进行独立预后分析（图7）：风险评分是TC 患者生存预后的独立因素（HR=8.882，95%CI：1.561～50.547，P=0.014），可独立于其他临床特征评估TC 患者生存预后。年龄（HR=1.078，95%CI：1.016～1.145，P=0.013）与TC 患者生存预后相关。

图7 TC 患者临床特征单因素和多因素Cox 回归分析森林图Fig.7 Forest plot of univariate and multivariate Cox regression analysis on clinical characteristics of TC patients

2.6 TC 差异PFRGs Nomogram 图构建及验证结果根据多因素Cox 回归分析结果，以年龄和风险评分作为因子构建Nomogram 图，预测TC 患者1、3 和5 年的OS 的概率。综合Nomogram 图中年龄及风险评分的总体得分，预测TC 患者的生存概率（图 8A），该 Nomogram 图 的 C-index=0.938（0.923～0.952），具 有 较 好 的 预 测 能 力。Calibration 分析结果显示：Nomogram 图预测TC 患者1、3 和5 年生存概率与实际观察结果吻合良好，见图8 B～D。

图8 TC 差异PFRGs 的Nomogram 图 和Calibration 图Fig.8 Nomogram and Calibration charts of TC differential PFRGs

2.7 TC 差异PFRGs 与蛋白共表达相关性分析结果Spearman 相关性分析基因共表达结果显示：有9 个编码基因与CD44 表达呈正相关关系（P<0.001），有95 个编码基因与ANXA1 表达呈正相关关系（P<0.001），有17 个编码基因与NR4A1 表达呈正相关关系（P<0.001）。分别绘制CD44、ANXA1 和NR4A1 相关系数前5 位的共表达热图（图9）：与CD44 最相关的前5 位基因依次是SFTA3 （cor=0.773）、SSBP2 （cor=0.770）、RNF146 （cor=0.755）、IKZF2 （cor=0.732）和ESYT3（cor=0.718）；与ANXA1 最相关的前5 位基因依次是DUSP6（cor=0.815）、SDC4（cor=0.812）、TMEM43 （cor=0.803）、BID （cor=0.793）和RAD23B （cor=0.783）；与NR4A1 最相关的前5 位基因依次是FOSB （cor=0.889）、NR4A3 （cor=0.845）、NR4A2 （cor=0.840）、CSRNP1（cor=0.837）和ZFP36（cor=0.826）。

图9 TC 差异PFRGs 共表达热图Fig.9 Heat maps of co-expression of TC differential PFRGs

PPI 结果显示（图10A）：PPI 网络中共有71 个节点，132 条边，即71 种蛋白质之间相互作用，共有132 组对应关系。拓扑分析结果显示：JUN、FOS 和ATF3 等 关 键 蛋 白 与CD44、ANXA1 和NR4A1 调控TC 铁死亡有关（图10B～D 和表2）。

表2 TC 差异PFRGs 关键蛋白拓扑分析Tab.2 Topological analysis on TC differential PFRGs key proteins

图10 TC 差异PFRGs PPI 网络图及拓扑分析图Fig.10 PPI network and topology analysis diagrams of TC differential PFRGs

2.8 TC 差异PFRGs GO 和KEGG 富集分析结果GO 和KEGG 富集分析结果显示：在满足AdjustedP<0.05 且q<0.2 条件下，TC 铁死亡主要富集于58 个BP、4 个MF 和5 个KEGG，主要与细胞周期的调控有关。聚焦于细胞周期BP、MF和KEGG（图11A）结果显示：TC 铁死亡主要富集于MAPK 活性失活、内肽酶活性的正向调节、凋亡通路的正向调节、ERK1 和ERK2 级联等生物进程；富集于MAPK 活性、p-MAPK 活性、生长因子结合和转化生长因子β 结合分子功能；主要富集于MAPK信号通路、Th17细胞分化和细胞凋亡等信号通路。

构建KEGG-靶点网络图（图11 B）：TC 铁死亡主要与MAPK 信号通路中双特异性磷酸酶（dual specificity phosphatase，DUSP1）1、DUSP5、DUSP6、Fos 原 癌 基 因 （Fos proto-oncogene，FOS）、白细胞介素1 受体辅助蛋白（interleukin 1 receptor accessory protein，IL1RAP）、Jun 原癌基因（Jun proto-oncogene，JUN）、MET 原癌基因（MET proto-oncogene，MET）、Ras 蛋白特异性鸟嘌呤核苷酸释放因子 1 （Ras protein specific guanine nucleotide releasing factor 1，RASGRF1）、转化生长因子α（transforming growth factor alpha，TGFA）、转 化 生 长 因 子β 受 体1 （transforming growth factor beta receptor 1，TGFBR1）和TNF受体超家族成员 1A （TNF receptor superfamily member 1A，TNFRSF1A）分子靶点的调控有关。

图11 TC 铁死亡GO 和KEGG 气泡图（A）及KEGG-靶点网络图（B）Fig.11 TC ferroptosis GO and KEGG bubble diagram（A）and KEGG-target network diagram（B）

3 讨 论

TC 是好发于中老年人，以高发病率和低死亡率为特点的内分泌恶性肿瘤［22］。目前，TC 根治术后肿瘤复发和不可切除患者缺乏有效全身治疗手段是TC 患者生存预后不良的主要原因［23］。为进一步提升TC 患者生存获益，突破临床治疗瓶颈，在肿瘤分子生物学研究的基础上，本研究结果显示：分子靶向药物可为碘难治性TC、不可切除性TC 和手术切除后复发TC 患者带来生存获益。随着肿瘤分子生物学研究及分子靶向药物研发的不断深入，索拉非尼、乐伐替尼和仑伐替尼相继问世，在新辅助治疗、靶向联合免疫治疗及姑息治疗方面极大地提高了TC 患者的生存获益［24］。铁死亡是由铁代谢、脂质代谢和氨基酸代谢共同介导的细胞程序性死亡模式，贯穿TC 发生发展的始终。研究［25］显示：针对肿瘤细胞铁死亡的分子靶向药物（索拉非尼）可明显延长肝癌、肾癌和食管癌等多种恶性肿瘤患者的生存期。本研究筛选出TC 差异PFRGs，分析其与TC 患者生存预后的关系，并构建预后风险评估系统及预测模型。

本研究通过TCGA、FerrDb 和GeneCards 数据库，差异分析鉴定出3 317 个在TC 中上调基因和3 456 个在TC 中下调基因，单因素Cox 回归分析鉴定出343 个差异基因与TC 患者的生存预后相关，其中包括3 个FRGs（CD44、ANXA1 和NR4A1）。结合基因表达谱及临床病理验证了CD44、ANXA1在TC 组织中明显上调，NR4A1 在TC 组织中明显下调，与差异分析结果相互佐证。CD44 是一种细胞黏附糖蛋白，在肿瘤进展和转移中发挥作用，被广泛用于TC 和其他恶性肿瘤的肿瘤干细胞标志物［26-27］，LI 等［28］采用双荧光素酶报告基因及免疫共沉淀实验证实miR-205-5p 能靶向GGCT，介导TC 细胞中CD44 的表达下调，抑制TC 细胞增殖、迁移和侵袭。ANXA1 是一种磷脂结合的膜定位蛋白，在TC 组织和细胞（SW579 细胞）中表达上调，可通过外泌体介导TC 和甲状腺滤泡上皮细胞（Nthy-ori3-1 细胞）之间的细胞通讯，促进TC 细胞及甲状腺滤泡上皮细胞的增殖、侵袭、上皮间质转化和肿瘤生长［29］。NR4A1 是核受体亚家族4A1组成员转录因子，入核调控下游靶基因转录，参与肿瘤细胞增殖、细胞凋亡和铁死亡等多个进程。JIANG 等［30］发 现：NR4A1 直 接 与 LEF1 的 启 动子区域结合，并导致与组蛋白乙酰化和 DNA 去甲基化的串扰，从而在转录上上调LEF1 表达，促进PTC 中下游生长相关基因的表达和肿瘤细胞增殖。

本研究采用Kaplan-Meier 曲线和time-ROC 曲线分析结果显示：CD44、ANXA1 和NR4A1 的表达与生存预后相关（P=0.048、0.005 和0.036），对TC 患者均具有良好的1、3 和5 年生存预测作用；采用单因素和多因素Cox 回归分析构建3 基因风险评分系统，计算TC 患者风险评分，结果显示：风险评分是独立于TC 其他临床特征的预后因子（HR=8.882，95%CI：1.561～50.547，P=0.014），风险评分越高，TC 患者1、3 和5 年生存预后越差（AUC=0.761、0.767 和0.722）；联合TC 患者临床特征构建预测TC 患者1、3 和5 年的生存概率的Nomogram 图（C-index=0.938），对TC 患者的生存获益具有较好的预测作用。

本研究中Spearman 相关性分析结果显示：CD44 与SFTA3、SSBP2、RNF146、IKZF2 和ESYT3 最 相 关；ANXA1 与 DUSP6、SDC4、TMEM43、BID 和RAD23B 最 相 关；NR4A1 与FOSB、NR4A3、NR4A2、CSRNP1 和ZFP36 最相关。联合3 个TC PFRGs 进行PPI 及拓扑分析结果显示：TC 铁死亡与JUN、FOS 和ATF3 等关键蛋白的相互作用有关。为研究TC 铁死亡介导的下游作用机制，本研究进行GO 和KEGG 富集分析，发现TC 细胞铁死亡主要与其共表达基因（DUSP1、DUSP5、DUSP6、FOS、IL1RAP、JUN、MET、RASGRF1、TGFA、TGFBR1 和TNFRSF1A）介导MAPK 信号通路，影响MAPK活性和p-MAPK 活性等分子功能，调控MAPK 失活有关。MAPK 信号通路与活性氧（reactive oxygen species，ROS）诱导的细胞铁死亡有关，细胞内过量的ROS 以脂质过氧化的形式介导铁死亡［31］，参与多种肿瘤细胞的增殖、迁移和侵袭过程［32-33］。

综上所述，本研究筛选得到的差异PFRGs（CD44、ANXA1 和NR4A1）与TC 发生发展和预后相关；成功构建了3 个基因联合的风险评分系统及预测模型，结合TC 患者临床特征能有效预测TC 患者1、3 和5 年生存期，TC 铁死亡的作用机制可能与介导共表达基因、激活MAPK 信号通路和调控相关分子功能有关。