基于整合生物信息学构建的肺腺癌内质网应激相关预后模型

熊雅俊，许爱国

(郑州大学第一附属医院 呼吸与重症三科，河南 郑州 450052)

原发性肺癌是我国最常见的恶性肿瘤[1]。非小细胞肺癌(non-small cell lung cancer，NSCLC)占新发肺癌的80%～85%[1]，肺腺癌(lung adenocarcinoma，LUAD)是NSCLC中最常见的病理类型，在年轻女性和从不吸烟者中的发病率不断升高。尽管以分子遗传学为基础的联合疗法改善了LUAD患者的预后，但由于缺乏早期识别标志物，61%的肺癌患者确诊时已进展为Ⅲ期或Ⅳ期[2]。据统计，NSCLC的5 a生存率只有23%，伴晚期转移性病变的LUAD患者的5 a生存率低于4%[2]。随着二代测序技术的跨越式发展，生物信息学的发展促进了基因表达谱技术在鉴定新型生物标志物方面的广泛应用[3]。内质网在蛋白质折叠、翻译后修饰、钙稳态和脂质合成中发挥重要作用[4]。多种外源性和内源性因素干扰内质网稳态，当错误折叠蛋白超过一定水平时，会触发内质网应激反应[5]。现有研究证明，内质网应激是许多癌症的标志。内质网应激的关键调节基因XBP1和GRP78在肺癌组织的mRNA和蛋白质中高表达，可促进肿瘤的侵袭和转移[6-7]。因此，内质网应激相关基因的异常表达可能对LUAD有预后价值，可作为潜在的治疗靶点。然而，影响LUAD进展和预后的内质网应激相关基因尚未完全鉴定。本研究基于癌症基因组图谱(The Cancer Genome Atlas，TCGA)数据库确定了190个内质网应激相关差异表达基因(differentially expressed gene，DEGs)，通过生物信息学分析筛选核心基因，建立内质网应激风险模型，以预测LUAD的预后。

1 资料与方法

1.1 数据下载及处理从TCGA网站(https://portal.gdc.cancer.gov/projects/tcga)下载了526例LUAD及59例正常肺组织的mRNA测序信息和临床数据，具有完整信息的有428例。GeneCards是1个整合了几乎所有已知人类基因信息的综合数据库，登录GeneCards网站(https://www.genecards.org/)，在搜索框输入“endoplasmic reticulum stress”，选择检索结果与检索关键词相关性≥7的833个基因。

1.2 筛选DEGs为了确定与LUAD预后有关的基因，使用R软件DESeq2包筛选患者和正常人组织间的DEGs。将标准设置为校正后的P<0.05，且差异倍数的绝对值≥1，利用R软件tinyarray包画热图和火山图，将结果可视化。最后，将DEGs与内质网应激相关基因取交集，即为LUAD差异表达的内质网应激相关基因。

1.3 富集分析为了阐明和DEGs相关的潜在的基因功能注释和富集通路，利用R软件clusterprofiler包进行基因本体(Gene Ontology，GO)和京都基因与基因组百科全书(Kyoto Encyclopedia of Genes and Genomes，KEGG)富集分析[8]，其中GO分析包括分子生物学功能(molecular function，MF)、生物学过程 (biological process，BP) 和细胞学组分三大类[9]。

1.4 Lasso和Cox构建模型利用R软件survival包对DEGs进行单因素Cox回归分析以确定和LUAD预后有关的内质网应激相关基因，再依次进行Lasso回归分析和多因素Cox回归分析进一步缩小目标基因范围。基于多因素Cox回归分析得到的回归系数和所选内质网应激相关基因的表达水平构建模型。

1.5 模型可视化绘制风险森林图可视化多因素Cox回归模型，为了评价模型的预测能力，绘制受试者工作特征(receiver operating characteristic，ROC)曲线预测模型的准确性。基于多因素回归分析的结果，将内质网应激风险特征联合临床指标绘制诺莫图，预测LUAD结局事件的发生。利用校正曲线评估诺莫图预测LUAD患者3、5 a生存率的准确性。校正曲线中的预测曲线和实际曲线的吻合程度越高，诺莫图的预测能力越强。

1.6 统计学方法采用R 4.2.0软件处理数据。采用Wilcoxon检验进行组间差异分析。采用单因素Cox回归模型计算风险比(hazard ratio，HR)和95%置信区间。多元Cox回归模型可提示内质网应激风险特征作为患者独立预后标志物的价值。所有P值均为双侧，P<0.05为差异有统计学意义。

2 结果

2.1 鉴定DEGs结果显示肿瘤组织中13 949个基因的表达与正常组织相比有显著差异，其中有3 262个下调基因，10 687个上调基因。火山图显示了DEGs在LUAD样本和正常样本之间的分布(图1A)。热图横坐标表示基因信息，纵坐标表示分组信息，图1B显示了正常组织和LUAD组织之间的基因组表达差异。使用主成分分析方法对样本进行降维处理，以探讨DEGs是否能区分LUAD样本和正常样本，结果显示正常样本和LUAD组织中DEGs的表达存在差异(图1C)。

A为DEGs的火山图；B为正常样本和肿瘤样本DEGs的等级聚类；C为DEGs的主成分分析。图1 正常组织与LUAD组织基因组图谱的差异

2.2 DEGs的富集分析为了确定DEGs的生物学功能，进行GO和KEGG通路富集分析。与生物学过程有关的DEGs显著富集于体液免疫反应、细胞识别及补体激活等，与细胞学组分有关的DEGs显著富集于质膜外侧面、含胶原的细胞外基质、免疫球蛋白复合物等，与分子功能相关的DEGs主要富集于信号受体激活剂活性、受体配体活性和门控通道活性等(图2A)。KEGG分析发现DEGs主要富集于脂质与动脉粥样硬化、内质网蛋白质加工和cAMP信号通路等(图2B)。

A为GO富集分析；B为KEGG通路富集分析。图2 DEGs的功能富集分析

2.3 内质网应激风险模型的构建和评估将上述DEGs和833个内质网应激相关基因取交集，得到190个DEGs。为了探索内质网应激在LUAD中的预后价值，构建风险评分模型以评估LUAD的内质网应激状态。首先，对190个基因进行批量单因素Cox回归分析，筛选出88个与LUAD预后显著相关的基因(P<0.01)。随后，利用Lasso回归算法筛选出回归系数非零的最有价值的预测基因(图3A和B)。绘制箱线图(图3C)和ROC曲线(图3D)评估两个λ值构建的模型的预测作用，结果提示模型可用，选择按照1 min建模筛选出的24个基因。为了进一步缩小目标基因的范围，利用多因素Cox回归分析方法构建最优模型，最终筛选出8个基因(图3E)。其中，MBTPS2、SEC61G、FURIN和PKP2是LUAD的危险因素(HR>1)，EIF2AK3、CAV3、SELENOK和NLRP1是保护因素(HR<1)。绘制ROC曲线评估内质网应激风险模型预测预后的有效性，如图3F所示，模型预测1、2、3 a总生存期(overall survival，OS)的曲线下面积(area under curve，AUC)分别为0.74(0.66～0.82)、0.75(0.69～0.81)和0.76(0.70～0.82)。

A为88个基因的最小绝对收缩和Lasso系数分布；B为在Lasso模型中通过10次交叉验证调整参数选择；C为预测结果的可视化；D为模型的预后预测能力评估；E为8个基因的风险森林图；F为模型预测LUAD患者1、2和3 a生存率的准确性。图3 预后模型的构建

2.4 内质网应激风险特征的功能注释根据MaxStat包确定的临界值0.92(图4A)，删除临床信息缺失者，将428例患者分为高危组和低危组。低危组患者表现出显著的生存获益(HR=0.28，95% CI：0.19～0.40)，见图4B。高危组中位生存时间为32.5个月，低危组中位生存时间为89.4个月。PKP2、MBTPS2、SEC61G和FURIN在高危组患者中的表达高于低危组，SELENOK、CAV3、NLRP1和EIF2AK3正相反(图4C)。为了探索内质网应激与LUAD相关的潜在机制，使用Hallmark通路基因特征分析方法鉴定两组富集到的生物学通路。结果显示高危组患者显著富集于mTORC1信号通路、G2M检查点和P13K/AKT/mTOR信号通路等和肿瘤发生发展相关的通路(图4D)。

A为将患者分为高危组和低危组；B为两组的Kaplan-Meier生存曲线；C为8个基因在两组的表达差异；D为两组的基因集变异分析。图4 构建riskscore特征

2.5 诺莫图的构建与验证多因素Cox回归分析显示，内质网应激风险特征和TNM分期与OS显著相关(图5A)。ROC曲线提示模型预测预后的有效性，如图5B所示，预测1、2、3 a OS的AUC分别为0.75(0.67～0.83)、0.75(0.70～0.81)和0.77(0.72～0.84)。利用TCGA数据集构建的内质网应激风险特征是LUAD的独立预后因素。为了在临床上更准确预测患者的病死率，整合内质网应激风险特征、年龄、性别和TNM分期，构建预测LUAD患者1 a和2 a OS的诺莫图(图5C)。校准曲线表明，患者3 a和5 a的预测生存率和实际生存率之间具有显著的一致性(图5D)。

A为多因素Cox回归模型；B为模型预测1、2、3 a生存率的能力；C为诺莫图的构建；D为诺模图的校准曲线。图5 构建诺莫图

3 讨论

肺癌是全球发病率和病死率最高的侵袭性肿瘤之一，NSCLC是最常见的肺癌类型，根据基因和分子特征又分为肺鳞癌和LUAD。近年来，诊断和治疗的进步极大地改善了LUAD的预后，但是由于肿瘤转移和复发，患者预后仍较差[10]。识别新的预后生物标志物以及建立预后模型具有重要意义。下一代测序技术的发展为解读LUAD的关键基因和表观遗传学的改变提供了机遇。

内质网在大多数蛋白质的合成和正确折叠中起关键作用[11]。失去调控的蛋白沉积导致错误折叠或未折叠蛋白的积累，引起内质网应激[12]。内质网应激是细胞在不利环境中生存的一种适应性机制[13]。癌细胞由于营养缺乏、缺氧和氧化应激等因素，容易诱发内质网应激[14]。内质网应激的持续激活使肿瘤细胞具有更强的致瘤、转移和耐药能力[14]。内质网应激参与肿瘤的发生和进展，了解内质网应激相关基因在LUAD中的作用，可能对基于内质网应激相关基因的临床结局预测和治疗靶点的开发至关重要。

既往文献曾发现多个内质网应激相关基因与肺癌发生和预后的关系。未折叠蛋白反应是内质网应激后的一系列适应性机制[15]。GRP78基因编码的蛋白是主要的未折叠蛋白反应调节因子，研究发现其在侵袭性、转移性和化疗耐药的肺癌中过表达[16]。GRP78单倍剂量不足能抑制肿瘤进展，延长患者生存期[17]。GADD34基因编码一种磷酸酶，负责内质网应激后恢复正常的蛋白质合成功能[18]。Lei等[19]通过在A549细胞系敲除GADD34，发现肿瘤细胞活性和增殖能力降低，肿瘤生长受到抑制。XBP1和AFT6基因编码的蛋白是内质网应激的关键调节分子，XBP1的mRNA和蛋白水平在肿瘤组织中过表达，与肿瘤侵袭性增加及上皮间质转化有关[6]。AFT6诱导表皮生长因子上调，刺激肿瘤血管生成，导致化疗后肿瘤复发[20]。本文旨在挖掘更多影响LUAD预后的内质网应激相关基因，有助于未来LUAD治疗靶点的开发，改善患者预后。

本研究从GeneCards网站检索并下载833个内质网应激相关基因。通过生物信息学分析，最终筛选出8个与LUAD预后相关的基因并构建预测模型。

PKP2基因编码一种结构蛋白，在调节蛋白结合、细胞连接和信号转导等方面发挥重要作用[21]。PKP2的高表达与卵巢癌[22]和LUAD[23]的增殖和侵袭有关，提示PKP2促进肿瘤进展，和本文得到的结果一致。表皮生长因子受体与肿瘤的发生发展有关，PKP2通过与表皮生长因子受体相互作用，促进其介导的信号通路活化，即细胞增殖、血管生成和侵袭能力的增加[24]。另有研究发现PKP2参与β-catenin介导的信号通路，是上皮间质转化的标志，PKP2可减少细胞间黏附，增强肿瘤细胞运动能力，这可能与其促进肿瘤迁移有关[25]。

NLRP1基因编码的蛋白通过形成炎症小体复合体，对诱导宿主防御过程中的炎症反应具有重要意义[26]，而炎症反应与癌症密切相关。Williams等[27]发现NLRP1在结肠癌患者中表达下调，在小鼠结肠肿瘤模型中，NLRP1可降低炎症驱动的结肠肿瘤的发病率和病死率。NLRP1表达下降与LUAD免疫细胞浸润程度降低及患者预后不良相关[28]，和本研究结果一致。NLRP1是一种抑癌基因，但其功能增益突变导致的异常NLRP1炎症小体激活与癌症的发生有关[29]。Zhai等[30]发现NLRP1通过促进炎症小体活化和抑制转移性黑色素瘤细胞凋亡促进肿瘤生长。Wei等[31]在乳腺癌MCF-7细胞中过表达NLRP1并建立了裸鼠移植瘤模型，发现NLPR1促进了乳腺癌的迁移、侵袭和生长。Zhong等[32]研究发现，NLRP1的胚系突变导致多发性自愈性掌跖癌和家族性慢性苔藓样角化症。

EIF2AK3基因编码一种内质网Ⅰ型跨膜蛋白，在恢复内质网稳态中发挥重要作用，本研究结果可知EIF2AK3的表达与LUAD预后良好呈正相关。Fei等[33]在LUAD细胞系中同样发现EIF2AK3基因水平降低。研究还发现在胰腺癌细胞中激活EIF2AK3/eIF2/ATF4信号通路能阻止肿瘤进展[34]。MBTPS2可作为肺癌患者的独立预后指标。Zhang等[35]发现MBTPS2在LUAD中高表达，与本研究结果一致。SEC61G基因编码的蛋白是SEC61易位复合体的1个亚基，在蛋白质折叠、修饰、易位和激活未折叠蛋白反应中起重要作用[36]。研究发现在A549细胞系中敲减SEC61G可抑制细胞增殖、迁移和侵袭能力，促进细胞凋亡[37]。SEC61G基因在肺癌中表达上调，与LUAD的不良预后显著相关[38]，与本研究结果一致。

本研究结果显示FURIN基因表达水平与LUAD预后良好呈负相关，CAV3和SELENOK基因正相反。FURIN基因编码的蛋白是钙依赖性前蛋白转化酶[39]。FURIN抑制剂通过下调迁移和凋亡相关蛋白的表达，抑制A549细胞的增殖和运动[40]。He等[41]研究发现FURIN在三阴性乳腺癌中高表达。CAV3基因在NSCLC中高表达，与患者预后不良相关[42]。敲减CAV3能抑制前列腺癌细胞的增殖、迁移和侵袭[43]。SELENOK编码一种内质网常驻蛋白，参与免疫细胞的抗氧化、钙通道调节和内质网相关降解等途径[44]。SELENOK基因表达降低与LUAD预后不良相关[45]。SELENOK的低表达与胃癌的增殖和侵袭能力增加有关，提示SELENOK基因抑制肿瘤进展[46]。

综上，LUAD高危组患者显著富集的通路与肿瘤的发生发展有关。基于多因素Cox回归分析的结果构建诺莫图，校正图提示诺莫图良好的预测潜能。因此，基于8个基因的内质网应激风险特征可以预测LUAD患者的OS，并有助于后续研究中选择最佳的治疗方案。然而，本研究也有一定局限性。首先，这8个基因在蛋白水平的表达和预后预测效果需要评估。其次，还需要进一步研究证实内质网应激风险特征在LUAD中的具体功能机制。此外，本研究结果需要结合湿实验进行验证。