基于网络铁死亡相关的长非编码RNA预测肝细胞癌患者预后研究*

布 凡，柴剑波，白 冰，王万宇，赵永厚

铁死亡相关(ferroptosis-related)是一种依赖于细胞内活性氧和细胞内铁的细胞死亡途径，不同于细胞凋亡和细胞自噬，与癌症发生、发展密切相关[1-4]。长链非编码RNA(long non-coding RNA，lncRNA)参与机体多个生理学过程，近期研究显示铁死亡相关lncRNAs与多种癌症患者的预后有关[5,6]。目前，铁死亡相关lncRNAs在HCC发病过程中的作用尚不清楚。本研究应用生物信息学方法获得影响HCC患者预后的铁死亡相关lncRNAs，为进一步探究HCC与铁死亡的关系以及评估HCC患者的预后提供新思路。

1 资料与方法

1.1 数据采集 自TCGA-HCC数据库收集50例正常人和374例肿瘤患者RNA测序数据。剔除临床资料不全的数据，实际纳入377例，其临床特征见表1。下载铁死亡数据库，包括铁死亡诱导基因、铁死亡抑制基因和铁死亡标志物，去除重复基因，共筛选出246个铁死亡相关基因。经Pearson相关性分析，计算铁死亡相关基因水平与lncRNAs的相关性，条件为P<0.001，相关系数绝对值>0.6，筛选出127个与铁死亡相关的lncRNAs。整理HCC患者临床和病理学资料，包括性别、年龄、肿瘤分级、临床分期、TMN分期、生存状况和生存时间，条件为错误发现率(false discovery rate，FDR)<0.01，差异倍数(fold change，FC)≥2时，进一步提取出73个差异水平的lncRNAs和80个差异水平的铁死亡相关基因。对差异基因进行GO/KEGG通路分析。

表1 TCGA数据库HCC患者临床信息

1.2 构建铁死亡相关lncRNAs预后风险模型 以P<0.01为筛选条件，采用单变量COX回归分析。根据铁死亡相关IncRNA(ferroptosis-related long non-coding RNA)水平和相应的风险回归系数，建立预后风险模型。计算每例患者基于该模型的风险值(risk score)，并以风险值的中位数将患者分为高与低风险组。计算公式如下:风险值=基因1)风险系数×基因(1)水平+基因(2)风险系数×基因(2)水平+…+基因(n)风险系数×基因(n)水平。

1.3 预测诺莫图(nomogram)分析 进行单样本基因集富集分析(single sample gene set enrichment analysis，ssGSEA)lncRNA参与的作用通路以预测预后。P<0.05设为有统计学意义，FDR-q<0.25为筛选条件，结合预后特征构建诺莫图，用于预测HCC患者1/3/5 a总生存率。

1.4 免疫与基因水平分析 应用CIBERSORT[7，8]、ESTIMATE[9]、MCP counter[10]、ssGSEA[11]和TIMER[12]5种算法评估铁死亡相关lncRNAs特征的高风险组与低风险组之间细胞成分和细胞免疫反应的异同。应用热图揭示不同算法下免疫应答的差异。

2 结果

2.1 差异水平的铁死亡相关lncRNA富集分析 对差异水平的铁死亡相关lncRNAs进行富集分析，绘制排名前30的通路(图1、图2)。GO功能富集分析得到762个条目，其中生物过程(BP)条目671个，细胞组成(CC)条目18个，分子功能(MF)条目73个。KEGG通路富集筛选得到56条信号通路。参与的BP主要包括细胞的化学应激反应、氧化应激反应、对营养水平的应答、对细胞外刺激反应；CC主要包括细胞顶端、粘着斑、细胞-底物连接、黑素体；MF主要包括氧化还原酶活性，作用于NADPH氧化酶复合物、有机阴离子跨膜转运蛋白活性、抗氧化活性、铁离子结合等。参与的相关通路主要有神经变性多种疾病通路、化学致癌作用-活性氧、癌症中的microRNA合成、血清素能突触、脂质和动脉粥样硬化、铁死亡、癌症的中心碳代谢等。

图1 差异水平的铁死亡相关基因KEGG富集分析

图2 差异水平的铁死亡相关基因GO富集分析

2.2 铁死亡相关lncRNAs预后分析 对73个差异水平的铁死亡相关lncRNAs进行单变量COX分析，确定了7个基于HCC差异水平的铁死亡相关lncRNAs(表2)。

表2 基于HCC差异水平的铁死亡相关lncRNAs

2.3 生存结果和多变量分析 Kaplan-Meier生存分析结果显示，高风险组总生存率显著低于低风险组(P<0.001，图3)。同时，ROC曲线显示预测3 a生存率的AUC为0.716，在预测HCC预后方面优于传统的临床和病理学指标(图4、图5)。通过风险分数(图6)和生存状态(图7)分析发现患者的风险评分与HCC患者的生存率成反比。风险热图显示，本研究发现的较多数lncRNAs与风险模型呈负相关，进一步排除了2个模型基因。差异IncRNAs预测1 a、2 a、3 a生存率的AUC为0.781、0.738和0.716(图8)。单变量COX分析显示，IncRNAs肿瘤分期(HR：1.68，95CI：1.37～2.06)，T分型(HR：1.66，95CI：1.37～2.02)和风险分数(HR：1.33，95CI：1.23～1.44)是HCC患者总体生存率的依赖性预后因素。对比多变量COX分析发现，IncRNAs肿瘤分期(HR：1.34，95CI：0.64～2.81)和T分型(HR：1.18，95CI：0.58～2.39)对HCC患者总体生存率的影响略有降低，而风险分数(HR：1.35，95CI：1.23～1.49)基本无变化。

图3 Kaplan-Meier生存分析

图4 风险预测的AUC

图5 DCA分析风险因素

图6 风险分数图

图7 生存状态图

图8 预测HCC患者1/2/3 a生存率

2.4 lncRNAs and mRNAs共水平网络和ssGSEA基因集富集分析应用Cytoscape软件可视化表现lncRNA和mRNA之间的共水平网络。低风险组与高风险组在年龄构成、性别、TNM分期方面无显著性差异，而低风险组肿瘤分期和分级多较高，高风险组则较低。结合临床和病理特征与铁死亡相关lncRNAs预后绘制诺莫图(图9)，预测1例55岁处于第3阶段的女性患者，其1 a、3 a和5 a病死率分别为0.368、0.681和0.852。高危组富集了多个肿瘤相关通路，如均质重组、癌症途径、mTOR信号通路和ERBB信号通路等。

2.5 免疫与差异mRNA水平分析基于CIBERSORT、ESTIMATE、MCP counter、ssGSEA和TIMER等5种算法绘制免疫反应热图。相关性分析表明，高风险组与低风险组细胞溶解活性、MHCI类分子、I型、II型INF反应存在显著性差异(图10)；免疫检查点显示，两组CD44、TNFRSF4和CD276等水平也不同(图11)。鉴于基于检查点抑制剂的免疫治疗的重要性，进一步分析发现两组患者之间YTHDF1、METTL3和METTL14等水平存在显著性差异(图12)。

图10 ssGSEA用于检测免疫细胞亚群和相关功能之间的关联

图11 高/低HCC风险人群免疫检查点水平比较

图12 高/低HCC风险人群m6A相关基因水平比较*P<0.05；**P<0.01；***P<0.001

3 讨论

在本研究，我们首先基于TCGA数据集确定了基于HCC差异水平的铁死亡相关lncRNAs特征，进一步探讨了肿瘤微环境下免疫浸润细胞和免疫检查点抑制剂在HCC患者预后中的作用。这项研究的发现揭示了铁死亡信号通路中潜在的生物标志物和治疗靶点，5个差异水平的铁死亡相关lncRNAs，即AC099850.3、LUCAT1、AL031985.3、MKLN1-AS和SREBF2-AS1，可能是HCC患者独立的预后因素。

LUCAT1位于细胞核和细胞质，可对诸多癌症起到调节功能，如肺癌[14]、肾透明细胞癌[15]、结肠癌[16]、肝细胞癌[17]等。对HCC患者进行回顾性分析，发现LUCAT1在HCC组织具有高水平，其水平活跃度通常与患者生存周期呈负相关。因此，将LUCAT1作为预测HCC患者预后的生物标志物。为了进一步确定LUCAT1对HCC的影响方式，在HCC异种移植模型的实验表明，通过上调LUCAT1水平，可以显著促进HCC细胞增殖、迁移和侵袭程度，在HCC的发展和转移过程中起到关键的作用。转录因子ETS原癌基因1(ETS1)在包括HCC在内的不同种类的癌症中发挥了致癌作用。假设MKLN1-AS通过结合miR-22-3p介导EST1水平，通过采集HCC患者癌组织，采用RT-qPCR检测MKLN1-AS水平，数据集显示MKLN1-AS主要位于HuH7和LM3细胞的细胞质中，HCC组织MKLN1-AS水平远高于正常组织。构建的人体肿瘤裸鼠异种移植模型表明，沉默MKLN1-AS能够抑制HuH7和LM3细胞增殖、血管生成、迁移和侵袭，证实了MKLN1-AS在HCC发病过程中作为致癌调节因子而发挥作用。基于多个数据库分析MKLN1-AS和Yes相关转录调节因子1(YAP1)在HCC组织水平，以确定MKLN1-AS在细胞中的定位以及MKLN1-AS对HCC患者的影响。结果表明，MKLN1-AS在HCC组织水平明显高于正常组织且主要位于HCC细胞质，MKLN1-AS高水平能够增强YAP1 mRNA的稳定性，进一步加速HCC细胞增殖、迁移和侵袭。在正常组织，MKLN1-AS可通过体内诱导YAP1水平而导致HCC的发生。因此，可以将MKLN1-AS作为YAP1的上游因子，用于HCC的诊断和预后判断。