肾细胞癌(RCC)是起源于肾小管上皮系统的恶性肿瘤，发病率居泌尿系统恶性肿瘤第二位，且并发症发生率每年呈上升趋势^[1]。肾透明细胞癌(Kidney renal clear cell carcinoma, KIRC)是最常见的肾细胞癌，其病因复杂多样，与遗传、吸烟、肥胖、高血压及抗高血压治疗等有关^[2]。KIRC早期症状不典型，大多数患者就诊时已是晚期，晚期患者放化疗疗效差，免疫治疗和免疫联合靶向治疗是发展的新趋势，然而仍存在诸多弊端，如费用昂贵、脱靶效应等^[3-4]。由于KIRC的病因复杂，肿瘤组织高度异质性等，患者的治疗效果仍不理想，因此急需寻找新的标志物，用于指导KIRC的治疗。

铁死亡是一种新的程序性细胞死亡方式，其特征在于铁依赖性的脂质过氧化物(LOS)积累使氧化还原反应平衡被破坏，最终导致细胞死亡^[5]。近年来，铁死亡与肿瘤之间的关系受到了广泛关注。研究表明GPX4是铁死亡发生过程中的关键调节器^[6]，一方面，通过抑制胱氨酸/谷氨酸逆转运蛋白体(System xc-)会引起谷胱甘肽(GSH)耗竭，进而使谷胱甘肽过氧化物酶4(Glutathione peroxidase 4，GPX4)失活；另一方面，铁死亡诱导剂可以直接抑制GPX4的活性。通过以上途径使GPX4失活，导致LOS的积累进而引发铁死亡的发生。越来越多的研究证实，肝细胞癌^[7]、胃癌^[8]、乳腺癌^[9]等对铁死亡敏感，但针对KIRC的相关研究较少。

由于铁死亡与肿瘤的关系受到许多铁死亡相关基因(Ferroptosis-related genes, FRGs)的调控，肿瘤组织中FRGs的表达可以用来评估患者的预后，例如，在肝细胞癌中，研成功构建相关预后模型，靶向FRGs可能是治疗肝细胞癌的一种选择^[10]。本研究通过分析TCGA数据库，获得KIRC的基因表达信息，进而分析样本中FRGs的差异表达，从而构建包含多个FRGs的模型去有效预测KIRC患者的生存情况，并通过分析风险评分模型与免疫状态之间的相关性，探讨潜在机制，为临床治疗提供诊疗依据以及寻找新的治疗靶点。

通过TCGA数据库(The Cancer Genome Atlas, https://portal.gdc.cancer.gov/)下载KIRC患者的mRNA表达谱数据，该数据集包括539例KIRC患者的肿瘤组织测序信息和72例正常组织测序信息，并利用R语言(R 3.6.3, https://www.r-project.org/)对数据进行提取和标准化处理。通过文献检索获得60个与铁死亡相关基因的信息^[10]，将表达矩阵中的基因信息与铁死亡相关基因取交集，获得新的铁死亡相关基因表达矩阵。由于TCGA数据库属于公开获取资源，本研究免于伦理委员会的批准。

使用R语言的“limma”包来识别肿瘤组织与正常组织中差异表达的FRGs，设置过滤条件：FDR(BH)矫正后的阈值P.adj＜0.05，对数差异表达倍数变化绝对值|log₂FC|＞1。

利用perl语言(perl 5.32.0, https://www.perl.org/)将铁死亡相关基因表达量与生存时间进行合并，删除临床信息不完整以及生存时间＜30 d的样本。利用“survival”包进行单变量Cox回归分析，筛选出具有预后价值的铁死亡相关基因；将具备预后价值的铁死亡相关基因同差异表达的铁死亡相关基因取交集，为防止过度拟合，利用lasso回归分析筛选掉相关性高的基因；之后进行多变量Cox回归分析，得到风险基因，构建预测生存时间的风险评分模型，并输出纳入模型基因的风险评分系数(Coef)。所构建的预后风险评分模型计算公式为：风险值(Risk score)=风险基因表达量1×coef 1+风险基因表达量2×coef 2+...+风险基因表达量n×coef n。

根据风险评分模型公式，计算513例患者的风险值，根据中位数值将患者分为低风险组和高风险组。运用Kaplan-Meier (K-M)法进行生存分析，明确风险值和生存时间之间的关系；采用时间依赖的偏移工作特征曲线(ROC)进行模型预测效果的评价；绘制高低风险组的风险曲线和生存状态图。利用单因素和多因素独立预后分析评估各个变量及风险值与预后之间的关系，判断风险模型能否作为独立预后指标。此外，利用“prcomp”包功能进行主成分分析(Principal Component Analysis, PCA)，还进行了t-SNE测验，将分组可视化，从而探索不同群体的分布情况。

根据铁死亡相关预后模型将患者分为高低风险两组，利用“limma”包对两组进行基因的差异表达分析，以此探究高低表达组之间生存差异的可能原因，设置过滤条件：FDR(BH)矫正后的阈值P.adj＜0.05，对数差异表达倍数变化绝对值|log₂FC|＞2。利用“clusterProfiler”“enrichplot”以及“ggplot2”三个包对差异基因进行GO(Gene Ontology)功能注释和KEGG(Kyoto encyclopedia of genes and genomes pathway, KEGG pathway)通路富集分析，设置P.adj＜0.05，Q.adj＜0.05，并绘制富集气泡图与条形图。

利用“gsva”包进行ssGSEA(Single sample gene set enrichment analysis, ssGSEA)分析16种免疫细胞浸润分数和13种免疫通路的活性。

对539例KIRC患者的肿瘤组织测序信息和72例正常组织测序信息进行铁死亡相关基因的差异分析，共得到21个差异基因(见图 1)，差异信息(见表 1)。

图 1(Figure 1)

图 1 差异表达的铁死亡相关基因 Figure 1 Differential expression of ferroptosis-related genes

表1(Table 1)

表 1 差异表达的铁死亡相关基因 Table 1 Differential expression of ferroptosis-related genes Gene conMean treatMean logFC P.Value FDR(BH) ACSL4 15 313.541 670 5 032.280 148 -1.605 524 1.18×10-36 6.50×10-35 AKR1C1 5 315.708 333 2 066.751 391 -1.362 897 2.55×10-09 4.53×10-09 AKR1C2 1 475.500 000 492.595 547 -1.582 728 8.78×10-05 1.21×10-04 ALOX5 708.000 000 2 809.239 332 1.988 358 1.57×10-22 5.75×10-22 ALOX12 73.305 556 156.350 649 1.092 791 7.59×10-18 2.20×10-17 ATP5MC3 19 642.333 333 9 230.107 607 -1.089 547 7.24×10-29 4.42×10-28 CBS 24.472 222 9.833 024 -1.315 438 1.38×10-16 3.78×10-16 CD44 4 693.055 556 11 508.252 000 1.294 069 7.07×10-14 1.49×10-13 CHAC1 617.208 333 124.230 055 -2.312 743 1.46×10-28 8.06×10-28 GLS2 68.402 778 33.345 083 -1.036 580 8.35×10-25 3.53×10-24 HMGCR 4 551.125 000 2 023.224 490 -1.169 566 1.54×10-30 2.50×10-29 MT1G 24 781.736 110 1 022.528 757 -4.599 064 8.52×10-30 6.69×10-29 PTGS2 2 080.791 667 786.341 372 -1.403 905 2.51×10-15 6.00×10-15 SLC7A11 71.833 333 223.278 293 1.636 118 5.90×10-21 2.03×10-20 FDFT1 8 185.180 556 3 540.884 972 -1.208 904 1.82×10-30 2.50×10-29 ACO1 12 368.138 890 5 520.211 503 -1.163 833 3.04×10-24 1.19×10-23 ACSL3 7 480.277 778 3 453.289 425 -1.115 121 5.58×10-29 3.84×10-28 PEBP1 79 704.555 560 30 176.905 380 -1.401 217 3.48×10-30 3.38×10-29 GOT1 13 892.986 110 5 944.525 046 -1.224 723 5.63×10-26 2.58×10-25 HMOX1 3 792.944 444 18 661.330 240 2.298 662 3.68×10-30 3.38×10-29 ACSF2 9 301.944 444 932.654 916 -3.318 117 1.92×10-32 5.28×10-31

表 1 差异表达的铁死亡相关基因 Table 1 Differential expression of ferroptosis-related genes

采用单因素Cox回归模型对铁死亡相关基因进行分析，获得28个与KIRC预后相关的基因(见图 2a)。将预后相关基因与差异的铁死亡相关基因取交集，共得到12个基因绘制Venn图(见图 2b)。根据单因素Cox回归分析结果，进行Lasso回归与多因素Cox回归分析，结果显示有10个基因被纳入模型(见图 2c、2d)。基于最佳λ值，构建出铁死亡相关基因的预后风险评分模型：风险值(Risk score)=(0.024 5)×ALOX5表达值+(0.126 0)×CBS表达值+(0.199 5)×CD44表达值+(0.218 3)×CHAC1表达值+(-0.295 9)×HMGCR表达值+(0.036 7)×MT1G表达值+(0.061 4)×SLC7A11表达值+(-0.080 7)×FDFT1表达值+(0.160 3)×PEBP1表达值+(-0.220 5)×GOT1表达值。

图 2(Figure 2)

图 2 铁死亡相关基因差异分析 Figure 2 Differential analysis of ferroptosis-related genes

ROC曲线分析(见图 3b)，显示1年曲线下面积(AUC)为0.739，2年为0.717，3年为0.713，均大于0.7，表明用上述模型预测KIRC患者预后的准确性高。K-M生存分析显示，高风险组与低风险组的生存有统计学意义(P=5.73×10^-13)，随着时间的推移，生存率逐年下降，其中高风险组生存率显著低于低风险组(见图 3c)。PCA和SNE分析显示，不同的风险组分布在两个方向上，表明10个基因的预后风险评分模型可以将高低风险组很好的区分开(见图 3d、3e)。

图 3(Figure 3)

图 3 预后风险评分模型评价 Figure 3 Evaluation of prognostic risk scoring model

计算患者的风险值，将其从低到高排序后依据中位值，分为高风险组和风险组。其中，高风险组有256例，低风险组有257例。生存状态图表明，高风险组死亡病例数比低风险组高，高风险组中总体生存率较短的患者多于低风险组，高风险组的ALOX5、CBS、CD44、CHAC1、HMGCR、MT1G、SLC7A11、FDFT1、PEBP1、GOT1表达量均高于低风险组(见图 3a)。

为评价预后风险评分模型是否能做为独立预后预测因子，进行单变量和多变量Cox回归分析，在单因素Cox回归分析显示年龄、grade分级、stage分期、T分期、M分期、N分期和风险评分与KIRC患者的OS显著相关(见图 4a)，其中风险评分的HR值为2.848，95%置信区间为2.124-3.819。在多因素Cox回归分析中，年龄、风险评分与KIRC患者的OS显著相关(见图 4b)，其中风险评分的HR值为2.173，95%置信区间为1.506-3.135。所以，风险评分可以作为评价KIRC患者的独立预后因素。

图 4(Figure 4)

图 4 单因素和多因素Cox回归分析结果 Figure 4 Results of univariate and multivariate Cox regression analyses

将差异表达的铁死亡相关基因进行富集分析，得到62个生物过程(BP)、5个细胞成分(CC)、4个分子功能(MF)。GO BP显示差异基因显著富集于经典途径的补体激活、免疫球蛋白介导的免疫应答、B细胞介导的免疫应答等免疫方面。GO CC显著富集于免疫球蛋白复合物、质膜、细胞外基质等方面。GO MF显著富集到抗原结合、免疫球蛋白结合等方面(见图 5a)。通过KEGG富集分析显示，差异表达的基因主要富集在IL-17信号通路上(见图 5b)。

图 5(Figure 5)

图 5 GO和KEGG分析的代表性结果 Figure 5 Representative results of GO and KEGG analyses

为进一步探索预后风险评分模型与免疫状态之间的相关性，使用ssGSEA量化了不同免疫细胞亚群、细胞功能、途径等，结果显示大部分免疫细胞的评分在高风险组和低风险组之间表现出显著差异(见图 6a、6b)，包括aDCs、巨噬细胞、Th1细胞、Th2细胞等。高风险组巨噬细胞、CD8+T细胞、APC共刺激等评分较高，而aDCs、Th1细胞、Th2细胞、T细胞共刺激等评分较低。表明高风险组抗原提呈过程较强，而细胞免疫、体液免疫等过程减弱。因此，高危患者抗肿瘤免疫减弱可能是其预后不良的一个原因。

图 6(Figure 6)

图 6 不同风险组间ssGSEA评分比较 Figure 6 Comparison of ssGSEA scores among different risk groups

铁死亡是一种全新的细胞程序性死亡，本质是铁离子依赖的脂质过氧化产物超量蓄积引起的以线粒体改变为主的氧化损伤，其存在于多种疾病的病理生理过程中，研究发现诱导细胞铁死亡可以杀灭肿瘤细胞进而抑制肿瘤的生长^[11]。

通过基因差异表达分析得到28个与KIRC患者预后相关的基因，根据单因素CoX回归分析、Lasso回归分析最终筛选得到10个与患者预后相关的基因(ALOX5、CBS、CD44、CHAC1、HMGCR、MT1G、SLC7A11、FDFT1、PEBP1、GOT1)，并以此为基准构建预后风险评分模型。

预后模型所用到的基因大致分为四类：脂质代谢(PEBP1、FDFT1、HMGCR、ALOX5)，抗氧化代谢(SLC7A11、CBS)，能量代谢(GOT1)和其他(CD44、CHAC1、MT1G)^[12]。ALOX5(花生四烯酸5-脂氧合酶)基因编码脂氧合酶，它是催化花生四烯酸转化为白三烯的关键酶，主要在白细胞中表达。研究表明该基因及其产物出现在包括乳腺癌、前列腺癌在内的多种癌症中，其代谢物在肿瘤细胞增殖、抑制凋亡以及转移中具有调节作用^[13]。标准培养条件下的前列腺癌细胞不断生成ALOX5的代谢产物，阻断这些代谢产物活性可导致癌细胞凋亡^[14]；而ALOX5基因多态性于肺癌风险增加相关^[15]。HMGCR(3羟基-3甲基戊二酰辅酶A还原酶)基因、FDFT1(维生素K2二磷酸法呢酰转移酶1)基因参与胆固醇的合成代谢，是调控脂代谢的重要基因。HMGCR在多种癌症中都有表达，在胶质瘤组织中通过激活YAP的转录活性增强胶质瘤细胞的成瘤性，而胃癌组织中则通过激活hedgehog/Gli1信号促进胃癌细胞的生长、迁移，下调该基因表达则可抑制癌细胞的生长、迁移甚至避免肿瘤的发生^[16]。CBS(胱硫醚-β合成酶)基因参与同型半胱氨酸转硫途径的第一步反应，将同型半胱氨酸(Hcy)不可逆转化为半胱氨酸和α-酮丁酸，其催化作用需要维生素B6参与。胱硫醚-β-合成酶活性降低将引起高同型半胱氨酸血症，进而导致多器官组织损伤。转硫途径不仅移除了体内过多的Hcy，而且对体内氧化还原平衡有重要的维持作用。研究表明，正常结肠黏膜和增生性结肠息肉中CBS表达相对较低，肿瘤组织中高表达的CBS通过H2S依赖性调节癌细胞生长、增殖和迁移^[17]。

近几年来，随着细胞铁死亡与肿瘤研究热点的兴起，免疫浸润相关分析作为预后指标的研究也逐渐进入人们的视野。基于不同风险组之间的差异，我们进行了富集分析，意外地发现许多免疫相关的生物过程和途径得到了丰富。推测铁死亡可能与肿瘤免疫有密切联系是合理的。为了进一步探索风险评分模型与免疫状态的相关性，我们量化了不同的免疫细胞亚群、细胞功能及途径等得出：高风险组巨噬细胞、CD8+T细胞、APC共刺激等比例较高，而aDCs、Th1细胞、Th2细胞、T细胞共刺激等比例较低。表明高风险组抗原提呈过程较强，而细胞免疫、体液免疫等过程减弱。较高的风险评分与抗肿瘤免疫损伤有关，因此，高危患者抗肿瘤免疫减弱可能是其预后不良的一个原因。

以肿瘤与细胞铁死亡间的相关性为前提，筛选到10个与患者预后显著相关的基因并以此构建风险评分模型，为KIRC患者的预后提供了新思路、新方法，具有一定的新颖性；本研究所得到的基因可以作为KIRC基础研究与治疗的可靠靶点，有一定的应用价值；不足之处在于，没有对所得到的模型进行验证，从而提高评分模型的可靠性；没有进行动物、细胞水平进行功能验证，缺乏大规模多中心的临床试验进行检验。