联合免疫分型构建肺腺癌铁死亡相关长链非编码RNA风险评分模型①

于海洋 潘跃银 （中国科学技术大学附属第一医院肿瘤化疗科，合肥 230001）

肺癌是我国最常见同时也是致死率最高的恶 性肿瘤，五年生存率仅为4%～17%，其中又以肺腺癌（lung adenocarcinoma，LUAD）所占比例最高，因此寻找LUAD 的治疗新靶点是当务之急［1］。铁死亡是一种新型的非凋亡细胞死亡模式，其主要机制为铁依赖性的脂质过氧化损伤，这种导致线粒体功能障碍和细胞内的毒性脂质过氧化的机制在抑制癌症生长和发展过程中起关键作用［2-3］。lncRNA 是长度大于200 个核苷酸，不具有或仅有少许蛋白编码能力的RNA 序列［4］。lncRNA 还参与各种生物调控过程，包括与肿瘤的发生、发展和转移有关的过程，如LINC00618 在人类白血病细胞中低表达，其通过提高脂质ROS 和铁水平，同时抑制SLC7A11 表达来促进铁死亡发生［5-6］。然而，铁死亡相关lncRNA 在LUAD 中的意义尚未明确，本研究依据LUAD 患者铁死亡与lncRNA 之间的相关性，结合免疫分型构建新型风险评分模型以评估LUAD患者的预后。

1 材料与方法

1.1 数据收集与处理 TCGA 数据库（https：//portal.gdc.cancer.gov/）下载LUAD 患者的转录组数据和临床信息，包括59 例非肿瘤组织和535 例肿瘤样本。Perl 软件对数据进行预处理以及lncRNA 表达信息和相关临床病理参数的提取。下载GEO 数据库中GSE31210 数据集中的表达数据及临床信息以供外部验证。采用R软件通过共表达分析铁死亡相关基因与lncRNA 之间的关系，以筛选出与铁死亡相关的lncRNA。结合临床数据，利用单因素Cox 回归分析法及Kaplan-Meier（K-M）生存分析法筛选出与LUAD预后有关的铁死亡相关lncRNA。

1.2 一致性聚类与免疫分型 基于得到的预后相关lncRNA 的综合表达信息，使用R 软件中的“ConensusClusterPlus”包对肿瘤样本的不同亚组进行了一致性聚类鉴定，用K-M 生存分析法绘制亚组间的生存曲线进行判定，并绘制与临床特征相关性的热图。运用CIBERSORT 软件（https：//cibersort.stanford.edu/）对所下载的LUAD 患者mRNA 数据处理后可得到肿瘤微环境内浸润的免疫细胞，通过转录组的表达数据可反映免疫细胞的浸润情况［6-7］；结合患者风险评估值与浸润免疫细胞进行Pearson 相关分析，找出与患者肿瘤浸润免疫细胞表达情况的相关性［8］。运用R软件中“estimate”包评估每个患者肿瘤样本免疫细胞和基质细胞的浸润水平（即免疫评分和基质评分）。

1.3 建立风险评估模型 筛选出预后信息不完整的31 例样本后，LUAD 肿瘤样本（n=504）被随机分为训练集（n=252）和验证集（n=252），随后使用单变量Cox 回归分析来检验训练集中铁死亡相关lncRNA 表达水平与总生存期（OS）之间的关系。筛选出符合条件的候选铁死亡相关lncRNA（P<0.01）。此后，对数据进行LASSO 回归，利用R 中的“GLIMT”包筛选意义最显著的预后因素，得出与预后有关的铁死亡相关lncRNA，并计算其风险评分。最终得到构成风险评分模型的评分公式如下所示：风险评分∑i N= 1（Ei*Ci），其中N 表示构建风险评分模型的铁死亡相关lncRNA 数目，Ci 表示铁死亡相关lncRNA 系数，Ei 表示铁死亡相关lncRNA 的表达水平。

1.4 模型评价与验证 依据风险评分模型计算训练集中每个患者的风险得分。随后依据训练集风险评分值中位数，将训练集与验证集患者各分为高风险组（风险值>中位数）和低风险组（风险值≤中位数）。使用R 软件绘制高低风险组K-M 生存曲线以对比两组患者生存差异，绘制ROC 曲线并用曲线下面积（AUC）来评估模型预测患者预后的准确性。然后分别对训练集和验证集患者的风险评分与其他临床因素（年龄、性别、分期）进行单因素及多因素Cox 回归验证，以考察其是否为独立预后因子。随后通过临床分组模型进行进一步验证。最后利用来自GEO数据库中的样本进行外部验证。

2 结果

2.1 预后有关铁死亡相关LncRNA 筛选 筛选出预后信息不完整的31例样本后，本研究汇总了来自TCGA 数据库中504 例肿瘤样本、59 例非肿瘤样本以及GEO 数据中226个LUAD 患者的临床基线资料（表1）。利用R 语言共表达分析FerrDb 数据库（http：//www.zhounan.org/ferrdb/）259 个铁死亡相关基因与lncRNA，以筛选出与铁死亡相关的lncRNA。同时结合LUAD 患者的临床信息，单因素Cox 回归分析最终筛选11 个与LUAD 预后有关的铁死亡相关lncRNA（表2）。

表1 患者的临床基线数据Tab.1 Clinical baseline data of patients

表2 单变量Cox 回归鉴定11个铁死亡相关lncRNA 的P 值和风险比Tab.2 P values and hazard ratios of 11 ferroptosis-related lncRNA identified by univariate Cox regression

2.2 一致性聚类与免疫分析 基于R 软件中的“ConensusClusterPlus”包对504 个肿瘤样本的不同亚组进行一致性聚类分析（图1A～C）。结果输出为K=2～9 个子组，当K=2 时为最优的划分方式。根据最稳定的K 值，所有患者被成功地分为两个亚组。其中，亚组1（Cluster1）代表低水平的基因表达，而亚组2（Cluster2）代表较高水平的基因表达。差异表达基因的总体生存分析表明，Cluster1 的生存期明显延长（P<0.001，图1D）。

图1 一致性聚类分析与鉴定Fig.1 Consistent cluster analysis and identification

肿瘤微环境差异分析显示（图2A～C），在Cluster1 中免疫细胞和基质细胞含量明显高于Cluster2，Cluster1 肿瘤纯度显著低于Cluster2，这一结果与上述生存分析结果也相互佐证。通过CiberSort 方法评估了一致性聚类中每个患者22 种不同免疫细胞的相对比例，并对这两个风险组产生的CiberSort 输出进行对比汇总（图2D），包括静息树突状细胞、静息肥大细胞、静息CD4+记忆T 细胞、幼稚B 细胞、中性粒细胞、CD8+T 细胞、辅助滤泡T 细胞在内的免疫细胞在不同的风险组中富集。结果表明静息树突状细胞、静息肥大细胞、静息CD4+记忆T 细胞在Cluster1 中显著高表达，幼稚B 细胞、中性粒细胞、CD8+T 细胞和辅助滤泡T 细胞在Cluster2 中明显高表达（图2E）。

图2 肿瘤微环境与免疫细胞浸润Fig.2 Tumor microenvironment and immune cell infiltration

2.3 构建风险评估模型 使用单变量Cox 回归分析来检验训练集中铁死亡相关lncRNA 表达水平与OS之间的关系，从而筛选出符合条件的候选铁死亡相关lncRNA（P<0.01）。此后，对高维数据进行LASSO 回归，利用R 软件中的“GLIMT”软件包选择意义最显著的预后因素。最终得出9个与预后有关的铁死亡相关lncRNA，分别为：Z97989.1、AL445524.1、AL391807.1、AL161431.1、AC024075.1、AC246787.2、FLG-AS1、LINC01352、AC090559.1，依据模型公式计算各lncRNA风险评分，风险评分值=EZ97989.1 ×（-0.940 467 387 287 318）+EAL445 524.1×0.000 879 959 066 730 395+EAL391807.1×（-0.657 561 513 468 988）+ EAL161431.1×0.00 416 410 759 392 408+EAC024 075.1×（-0.0 111 194 380 084 745）+EAC246787.2×（-0.137 501 955 859 741）+EFLGAS1×1.94901939831152+ELINC01352×（-1.83 848 358 454 151）+EAC090559.1×（-0.159 643 403 496 863），其中E表示铁死亡相关lncRNA表达水平。

2.4 风险评估模型的评价与验证 利用“survival”包对高低风险组患者的生存进行差异分析，并绘制风险评分曲线及生存状态关系图（图3A）。随着风险评分的升高，相较于低风险组而言，高风险组患者的病死人数明显增多。生存分析表明，与高风险组相比，低风险组患者OS 显著延长（图3C）。为了进一步评估此模型的特异度和灵敏度，通过“survivalROC”程序包来绘制模型的ROC 曲线。ROC 曲线显示：训练集中1 年和5 年AUC 值分别为0.708、0.758（图3E），表明该模型具有较好的预测效能。通过单因素以及多因素Cox回归分析探讨了风险评分模型是否是LUAD 患者预后的独立预测因子。在除其他临床病理特征（如肿瘤分期、年龄以及性别）的干扰后，此模型具有较好的预测效能。

图3 风险评分模型预后评估Fig.3 Risk scoring model prognostic assessment

依据同一风险评分值，将验证集患者分为高、低风险两组。与从训练集得出的结果一致，随着风险评分的升高，相较于低风险组，高风险组患者的预后明显更差（图3B）；且高风险评分患者组的中位OS 与低风险评分组比要更短（P=0.001，图3D）。验证集中1 年和5 年AUC 值为0.659、0.636（图3E、F），这表明本预测模型在验证集也具有良好的预测效能。对临床因素和风险评分值与预后关系的分析表明，与LUAD 患者的肿瘤分期、N 分期、免疫分组均与风险评分相关（P<0.05），且风险评分随着患者肿瘤淋巴结分期增加，这说明此风险评分预后模型与LUAD 淋巴结转移存在一定关联。R 软件“survminer”包根据临床分组从而对模型进行进一步验证，如图4 所示，以年龄、性别、分期及T、N、M 分期作为分组依据，本风险预测模型同样适用于临床分组后的模型。依据同一风险评分值，将来自GEO数据库226 例LUAD 患者分为高、低风险两组，结果表明高风险组患者的预后明显差于低风险组（P<0.001，图5A），1 年、3 年及5 年ROC 曲线AUC 值分别为0.61、0.69、0.77（图5B），表明本预测模型在外部验证集同样具有良好的预测效能。

图4 风险评分值与临床因素Fig.4 Risk score values and clinical factors

图5 GSE31210数据集外部验证Fig.5 External validation of GSE31210 dataset

对免疫细胞表达情况和免疫预后RiskScore 风险值进行Pearson 相关性检验，由图6 可知，巨噬细胞M0、静息肥大细胞、活化肥大细胞、静息树突状细胞和活化树突状细胞等肿瘤浸润免疫细胞与患者免疫预后RiskScore风险值密切相关（P<0.05），活化树突状细胞、巨噬细胞M0、活化肥大细胞、活化NK细胞和辅助滤泡T 细胞呈正相关（R>0）；记忆B 细胞、静息树突状细胞、巨噬细胞M2、静息肥大细胞、单核细胞和静息CD4+记忆T细胞呈负相关（R<0）。

图6 预后风险评分与免疫细胞浸润Fig.6 Prognostic risk score and immune cell infiltration

3 讨论

精准的疗效预测和预后判断方法是临床防治LUAD 的重中之重。目前以临床特征或血清分子标志物（CEA、CA125 等）作为预后判断的方式较为局限，高通量测序技术的广泛应用则为癌症诊断、治疗和预后研究提供了越来越多的测序数据［9-10］。lncRNA 指可通过影响邻近的基因或其他染色体上的远距离基因来达到调控编码基因表达目的并且长度超过200个核苷酸的非编码RNA 序列［1］。既往研究提示在LUAD 中有很多异常表达的lncRNA（如lncRNA MIR31HG、lncRNA H19、lncRNA HOXCAS3 等），作为LUAD 患者特异性的生物学标志物在诊断、预测预后等方面起一定作用［11-13］。CAI等［14］研究表明，在裸鼠肺癌肿瘤模型中，过表达转铁蛋白受体1（transferrin receptor 1，TFR1）可促进肺癌细胞对铁吸收，从而促进肿瘤细胞的增殖。已有研究表明通过构建铁死亡相关基因风险评分模型可以预测肺腺癌患者的预后［15］；此外，LUO 等［16］通过构建铁死亡相关lncRNA 预后模型可在一定程度上预测头颈部鳞状细胞癌患者的预后。然而，目前为止尚未有预测肺腺癌患者的铁死亡lncRNA 风险评分模型。因此，构建相关风险评分模型有助于LUAD 患者的预后评估。

本研究经过对TCGA 数据库的LUAD 样本的汇总和筛选，最终得到504例符合标准的样本数据，并将其随机分为训练集与验证集，最终得出由9 个铁死亡相关的lncRNA 构成的风险评分模型，结果表明，AL445524.1、AL161431.1 和FLG-AS1 的高表达均与患者的 OS 呈正相关；Z97989.1、AL391807.1、AC024075.1、AC246787.2、LINC01352和AC090559.1 高表达与较差的预后相关。先前已有研究发现AL161431.1 及LINC01352 分别与肝癌和子宫内膜癌的发生发展联系紧密，且通过构建lncRNA 预后模型可在一定程度上预测患者的预后［17-19］。无论是单因素分析还是多因素Cox 回归分析均表明该风险预后模型是影响LUAD 患者预后的独立危险因素（P<0.001）。该风险模型对训练集和验证集中1年和5年AUC 值分别为0.708、0.758；验证集中1 年和5 年AUC 值分别为0.659、0.636。通过一致性聚类使用多种不同的聚类方法，从而找到一种比每种单独的方法更合适的聚类方法对肿瘤进行分组，以便有助于为每个患者制定个性化的治疗方法。根据基因表达水平，使用R 软件的“ConensusClusterPlus”包 将 数 据 分 为Cluster1 与Cluster2。生存分析表明，Cluster1 的生存时间明显延长，提示生存时间与铁死亡相关基因的综合表达水平有关。

然而，本研究仍存在一定程度的不足，本研究是一项基于TCGA 数据库的回顾性研究，有可能产生一定范围内的偏差。此外，由于缺少合适的数据集，本研究尚未在其他数据库进行深入的验证；且本研究是一项基于高通量测序结果的回顾性研究，并未进行进一步的功能研究和体内外实验。因此需要对本风险预后模型进行更加深入的探索和后续研究。