基于生物信息学分析和实验验证筛选影响肝癌患者预后的铁死亡调控基因

姜宇朗，王 铮，孙明瑜

（1.上海中医药大学附属曙光医院，上海中医药大学肝病研究所肝肾疾病病证教育部重点实验室，上海 201203；2.上海中医药大学，上海 201203）

原发性肝癌是目前我国发病率第四和死亡率第二的恶性肿瘤，严重威胁着我国人民的生命和健康安全［1］。由于肝癌起病并无特异性的临床症状，大部分患者在确诊时往往已经处于中晚期，丧失了最佳的手术治疗机会［2］，导致五年生存率较差，因此早诊断早治疗是肝癌防治的一大目标，对于肝癌诊断预后具有良好敏感性和特异性的预后生物学标志物仍亟待开发［3］。

铁死亡是近年来发现的一种依赖于铁离子的调控且有别于凋亡、焦亡、程序性坏死、自噬的新型调控性细胞死亡［4］，其主要的细胞形态学变化特点是线粒体的改变、包括线粒体的皱缩、膜密度的增高、嵴的减少等，生化学方面改变主要包括脂质过氧化物的蓄积，活性氧（reactive oxygen species，ROS）的增多，最终对磷脂膜造成损伤并使细胞死亡［4，5］。铁死亡作为一种代谢性细胞死亡方式受到多种代谢途径的调控，其中包括氧化应激、铁离子代谢、线粒体能量代谢和胱氨酸代谢、多不饱和脂肪酸的合成以及葡萄糖的代谢密切相关［5，6］。研究表明铁死亡在肝癌的发生、发展中都扮演着重要的角色，诱导肝癌细胞铁死亡有望发展为肝癌治疗领域的一颗新星［7，8］。

研究发现，铁死亡相关基因在癌症的预后情况和预防诊断中也存在着重大的价值，有研究发现17个铁死亡相关的长链非编码RNA 可以用来预测结肠癌的转归［9］，8 个在胃腺癌患者中差异表达的预后FRGs 可以预测胃腺癌患者的生存预后状态［11］。但是诊断肝癌预后转归的差异铁死亡基因仍然未被有效发掘，本文旨在通过生物信息学和实验验证的方式，找出影响肝癌患者预后的差异FRGs，从而为临床肝癌预后新的生物标志物的开发提供线索，并为肝癌铁死亡疗法发掘新的治疗效果判断的标志物。本研究的流程和具体方法见图1。

图1 生物信息学分析的工作流程示意图Fig 1 Schematic diagram of the workflow of bioinformatics analysis

1 材料和方法

1.1 基因表达量的下载和铁死亡相关基因的表达量的整理

肝癌患者的基因测序结果从癌症基因组图谱计划（The Cancer Genome Atlas， TCGA）（https：//portal.gdc.cancer.gov/）数据库下载获得，结果共获得422 个数据集，其中包括374 例肝癌患者和50 例正常样本。其中铁死亡相关基因从Ferrdb 数据库（http：//www.zhounan.org/ferrdb/）获取得到。将FRGs 和TCGA 数据库基因表达相关联，得到FRGs 在每个样本中的表达数据。

使用limma 包将肝癌组织和正常组织的基因表达量进行非参数检验，设定差异条件为错误发现率（false discovery rate， FDR）＜0.05 并 且| log2Fold-Change|＞1，获得在肝癌样本和正常组织中差异表达的相关FRGs。

通过综合分析患者的生存时间和生存状况，并结合FRGs 的差异表达，使用survival 包进行单因素COX 回归分析，从而确定与肝癌患者预后有关的FRGs（P＜0.05）。

1.2 影响肝癌患者预后的铁死亡相关差异表达基因的确定和PPI 网络的构建

首先使用venn 包来筛选出差异表达的FRGs和预后生存相关的FRGs 的共同部分，这样可以获得影响肝癌患者预后的差异FRGs。然后，使用pheatmap 包绘制出这些差异表达的预后FRGs 的热图。最后，通过survival 包绘制相应的森林图。将上述基因导入蛋白互作网络（protein protein interation network， PPI）工具（http；//string-db.org/）可视化蛋白之间的互作关系。为进一步体现蛋白质之间的表达关联性，计算出核心蛋白的表达相关性并使用igraph 包绘制蛋白质互作网络图。

1.3 风险模型的构建及其验证

利用glmnet 包中的LASSO 回归比例风险模型，构建模型并获取基因风险系数，计算风险评分。根据风险值中位数，将样本分为高低风险组。通过survival 和survminer 包，比较两组生存和预后差异。使用timeROC 包绘制时间依赖ROC 曲线，评估模型预测准确性，利用pheatmap 包绘制风险曲线。

1.4 PCA 和tSNE 降 维 分 析

使用Rstne 包进行主成分（principal component analysis， PCA）和 T 分 布 随 机 邻 嵌 入 分 析（t-distributed stochastic neighbor embedding， tSNE）观察患者根据高低风险模型的分组情况的分布。

1.5 独立预后因素分析

将患者的肿瘤分期、肿瘤等级、年龄、性别、风险评分使用survival 包进行单因素和多因素COX 回归分析，独立预测因素分析，并绘制单因素和多因素对生存时间和生存状态的影响的森林图。

1.6 风险差异基因分析

使用limma 包设定错误条件发现率为FDR 值小于0.05 且| log2FoldChange|＞1，将风险值分别属于高低组的人群进行基因差异表达分析，得到各个样本在TCGA 数据库中高低风险组中的基因的平均表达量。

1.7 GO 和KEGG 富集分析和免疫相关分析

使用enrich-plot 包将风险差异基因进行GO 和KEGG 富集分析，借助ggplot 包使用气泡图将结果可视化。使用GESAbase 和GSVA 包进行免疫细胞和功能的免疫浸润相关分析，并用ggplot 包绘制高低风险组患者的箱线图。

1.8 差异的预后铁死亡基因的免疫组化验证

在人类蛋白质图谱数据库（https：//www.proteinatlas.org/）下载TOP20 基因在正常肝组织和肝癌组织的高分辨率免疫组化结果图。

1.9 差异的预后铁死亡基因的细胞层面验证

裂解正常肝细胞L02 和肝癌细胞HepG2，提取细胞总RNA，测定其浓度和纯度后，将RNA 逆转录为cDNA，以cDNA 为模板进行qPCR 扩增。以正常组作为对照，采用RQ 值来显示目的基因mRNA表达水平。比较正常肝细胞和肝癌细胞的目的基因的表达量。

1.10 差异的预后铁死亡基因组织层面验证

选取来自曙光医院的六名临床确诊为肝细胞癌的患者，手术切除样本组织后称取肝癌和癌旁组织各50 mg，匀浆机匀浆， 提取组织总RNA， 进行后续PCR 检测。比较正常肝组织和癌旁组织目的基因的表达量。

1.11 统计学方法

数据用GraphPad Prism 9 软件进行统计，两组间计量资料符合正态分布且方差齐时，选用用t检验，多组间比较用One-way-ANOVA 单因素方差分析，进一步多组间比较采用Tukey 检验。若数据不符合正态分布则采用非参数秩和检验。所有数据以（±s）表示，P＜0.05 为差异有统计学意义。

2 结果

2.1 差异的肝癌预后FRGs 的筛选结果

从TCGA 数据库共获得374 例肝癌样本和50例正常样本的RNA-seq 数据和对应的临床数据集。从ferrdb 数据库上共获得FRGs249 个，通过分析得到84 个差异表达的FRGs 和90 个预后相关的FRGs，对两者取交集可以得到42 个和肝癌预后相关的差异表达的FRGs（图2A）。依据上述42 个基因表达数据绘制热图（图2C）。同时筛选出P值排名靠前的20 个基因，其风险比率（hazard ratio，HR）及其95%置信区间如表1 所示，同时绘制森林图（图2B），确定肝癌的风险基因，尤其以ZFP69B、EIFS21、ATG3、ATG7等铁死亡相关基因的致病风险效应最为明显。

表1 单因素COX 回归分析得出的与肝癌预后相关的FGRs（P值排名前20）Tab 1 FGRs related to liver cancer prognosis obtained from univariate COX regression analysis（top 20 P value）

图2 筛选出的影响肝癌预后的FRGsFig 2 Screened FRGs affecting prognosis in liver cancer

2.2 蛋白质相互作用网络的构建结果

将42 个交集与肝癌预后相关的差异FRGs 导入STRING 数据库，剔除游离靶点，得到基于String数据库的蛋白质互作网络（图3A），其中相关性最高的核心靶蛋白主要有MAPK、NRAS、SRC、HRAS和CDKN2A。蛋白相关性作用网络（图3B）结果显示，ALB 与其他蛋白之间呈现出负向调控的关系，而其余基因则基本表现出正相关。

图3 蛋白互作网络的结果示意图Fig 3 Schematic representation of the results of the protein interaction network

2.3 风险预后模型的建立以及检验

LASSO 回归分析得到肝癌预后的风险模型，最 终 获 得11 个 基 因（G6PD、HRAS、SLC1A5、MT3、SRXN1、SLC7A11、ZFP69B、SLC2A1、STMN1、RRM2、SQSTM1）。根据风险值的中位数，将患者分为高、低风险组，结果显示高风险组患者的生存率和生存时间显著低于高风险组患者（P＜0.01）（图4A）。患者三年生存率的ROC 曲线面积均大于0.6，说明本研究建立的风险模型可信（图4B）。同时根据患者的风险值及其分组绘制患者人数和风险值的关系以及存活时间和患者人数的关系（图4C，D），可以看出，随着风险值的增加病人的死亡数明显增加。为了进一步研究高低风险组人群的分布情况，使用PCA（图4E）和t-SNE（图4F）可视化降维分析。最后再进行单因素（图4G）和多因素（图4H）COX 回归分析来检测能单独影响肝癌患者预后的因素，研究发现，肝癌分期和风险评分可以作为肝癌患者预后的独立预测因子（P＜0.01）。

图4 影响肝癌患者预后的模型的构建及效能检验和肝癌患者独立预后因素的分析Fig 4 Construction and efficacy test of the model affecting the prognosis of liver cancer patients with independent prognostic factors

2.4 风险差异基因的GO、KEGG 富集分析结果

使用limma 包获得差异FRGs 在高低风险组的患者的平均表达量，并做差异分析，最终得到风险差异FRGs，GO 分析表明风险预后FRGs 主要和体液免疫应答、细胞吞噬作用、淋巴细胞介导免疫等生物学过程有关，KEGG 富集分析显示风险差异基因参与包括细胞周期、PI3K/AKT 通路、脂质氧化、癌症中的蛋白多糖的表达等（图5A，B），聚类分析结果显示，这些肝癌的风险FRGs 和免疫的关系极为密切，这些基因可能在调控肿瘤的免疫微环境中起作用。

图5 差异预后FRGs 的Go 和KEGG 富集分析结果气泡图Fig 5 Bubble plots of Go and KEGG enrichment analysis results for differential prognostic FRGs

2.5 免疫浸润分析

进一步为了探究差异预后FRGs 和免疫功能的关系，通过单样本基因集富集分析（ssGSEA）确定样本高风险组和低风险组的免疫细胞浸润评分和免疫相关功能评分（图6）。相关评分显示，高低风险组患者的aCDs、Macrophages、NK 细胞表达差异具有统计学意义，在免疫功能方面主要和抗原呈递细胞的激活、免疫检查点、MHC 分子、副炎症、2 型干扰素的应答有关。

图6 高低风险组患者的免疫细胞和免疫功评分Fig 6 Immune cell and immunocompetence scores of patients in high and low risk groups

2.6 正常组织和肝癌组织的差异预后FRGs 免疫组化结果

从人类蛋白质图谱数据库中下载得到了16 个基因在肝组织和肝癌组织中的免疫组化结果（图7）（部分基因在数据库中并无蛋白质表达数据）。分析目的基因在两者之间的相对表达量发现G6PD、PRDX1、TXNRD1、ATG3 在肝癌组织中表达明显上 调，而SLC1A5、NRAS、EIF2S1、SLC38A1、ATG7在肝癌组织中表达明显下调，而其余基因表达量变化并无明显差异。

图7 TOP20 基因在正常肝组织和肝癌组织中的相对表达的免疫组化结果图Fig 7 Immunohistochemical results of TOP20 gene expression in normal and hepatocellular carcinoma tissues.

2.7 肝细胞与肝癌细胞的差异预后FRGs 的基因表达结果

为进一步验证上述生物信息分析结果的可靠性，在肝细胞和肝癌细胞中进行PCR 扩增反应后，和肝细胞比较肝癌细胞中部分差异FRGs 发生显著改变（图8），肝癌细胞中G6PD、STMN1、SLC7A1 1、TXNRD1表 达 水 平 显 著 增 加（P＜0.05），而SLC1A5、ATG3水平显著下调（P＜0.05）。

图8 基因在肝细胞和肝癌细胞差异预后FRGs 的相对基因表达量Fig 8 Relative gene expression of TOP20 gene in differentially prognostic FRGs of hepatocytes and hepatocellular carcinoma cells

2.8 肝癌和癌旁组织的差异预后FRGs 的相对表达结量

和癌旁组织比较，肝癌组织的G6PD、SLC7A1 1、TXNRD1水平显著增加（P＜0.05），而SLC2A1、SLC1A5、ATG3表达水平显著降低（P＜0.05）。这与细胞实验的结果基本一致（图9）。

图9 TOP20 基因在肝癌和癌旁组织差异预后FRGs 的相对表达量Fig 9 Relative expression of TOP20 gene in differentially prognostic FRGs of hepatocellular carcinoma and Paracarcinomatous tissues

3 讨论

基于生物信息学的分析结果发现G6PD、ZFP69B、SLC2A1、SLC1A5、HILPDA、MAFG、STMN1、SRXN1、SLC7A11、NRAS在肝癌预后的风险比当中占比较高。意味着这些基因有进一步开发为肝癌预后的标志物的价值。根据多个差异表达的FRGs 基因构建的风险预后模型可以较好的将患者区分为高风险和低风险人群，为肝癌的诊断预防提供参考，这些基因也可能是预测靶向铁死亡疗法治疗肝癌的有效性预测的潜在生物标志物［12］。进一步的结合免疫组化和细胞与组织的PCR 实验结果发现其中G6PD、SLC1A5、TXNRD1的变化趋势和生物信息学的推测结果一致。这提示在这些差异预后的FRGs 当中，这些靶标可能位于整个肝癌铁死亡调控网络的核心位置。

其中G6PD是葡萄糖-6-磷酸脱氢酶的编码基因，其主要功能是产生 NADPH，而NADPH 是抵御氧化剂和还原性生物合成反应中的关键电子供体［13］。G6PD 是葡萄糖代谢的关键酶，在肝癌患者中G6PD 往往处于高水平，这也意味着肝癌中葡萄糖高代谢状态，而葡萄糖高代谢又会诱导铁死亡的发生，G6PD也可以作为肝癌患者预后的独立风险因素来对患者进行定期监测［14］。因此G6PD有成为肝癌铁死亡预后基因的潜力，但仍然需要进一步基础实验和临床实验佐证。

SLC1A5作为人溶质载体家族的一员，是一种钠离子依赖性的中性氨基酸的转运体，可以转运包括和铁死亡密切相关的谷氨酰胺、甘氨酸、谷氨酸、异亮氨酸等［15］。肝癌细胞氨基酸高代谢状态是肿瘤高侵袭力的体现之一，代谢重编程是肿瘤细胞的特征性印记之一［16，17］。靶向SLC1A5可以阻止谷氨酰胺进入细胞，从而降低内肿瘤细胞内谷氨酰胺介导的高新陈代谢状态，在一定程度上“饿死”肿瘤以达到癌症治疗的目的［18，19］。

在经典的铁死亡防御体系中，胱氨酸通过细胞内唯一的胱氨酸/谷氨酸反向转运体进入细胞，随后在TXNRD1 酶的作用下迅速被还原为半胱氨酸，其作为谷胱甘肽的底物之一促进随后的谷胱甘肽合成［20，21］。最后GPX4利用谷胱甘肽作为底物来还原胞内过多的脂质过氧化物，从而抑制铁死亡的发生。因此TXNRD1可以间接反应胞内GSH 水平，也能预测患者对于铁死亡的敏感性，从而评估铁死亡药物的临床疗效［22］。

筛选出的三个基因都和肿瘤的代谢密切相关，尤其是氨基酸的吸收、分布、转化方面，而氨基酸代谢作为铁死亡三大组成部分之一和铁死亡的关系不言而喻。诸多氨基酸的代谢都参与了铁死亡的调控，而其中可能可以进一步发掘出诱导铁死亡治疗疾病的新靶点。但是本文并未进一步研究这三者之间内在的调控关系，它们之间是否存在某些内在或者外在的关联性需要进一步探究。

总而言之，本研究筛选出的和肝癌相关FRGs可能用来预测常规病理因素之外的肝癌预后，与铁死亡相关的这些特异的生物标志物也有进一步发展成为肝癌的诊断性指标的可能，同时针对这些靶点的铁死亡途径也可能为肝癌的治疗提供新的方向。鉴于铁死亡和免疫的关系密切，未来癌症治疗的研究重点应放在免疫阻断剂和靶向治疗结合诱导细胞铁死亡［23，24］。

作者贡献度说明：

姜宇朗设计文章的思路与框架、撰写论文；王铮参与生信分析；孙明瑜指导论文撰写并审校。

所有作者声明不存在利益冲突关系。