肝细胞癌免疫相关基因和lncRNA联合预后模型的构建及验证

王威，程倩倩，周雪丽，季文斌，吕振宇，杨燕

作者单位：蚌埠医学院第一附属医院肿瘤内科，安徽 蚌埠233004

肝细胞癌（hepatocellular carcinoma，HCC）是常见消化系统恶性肿瘤，全球每年约有91万新发病例和83万死亡病例，在所有癌症中居第6位及第3位［1-2］。尽管近年来HCC新药不断涌现，尤其是免疫治疗或免疫联合治疗较大地改善了晚期HCC的整体预后，但有效疗效预测和预后判断生物标志物的缺乏在一定程度上限制了人群的获益。传统临床分期虽是指导HCC治疗及判断预后的重要因素，但一些临床特征相同的病人，由于分子层面异质性可能预后不同［3］。近年来，大量证据显示包括HCC在内的肿瘤免疫微环境与病人的生存关系密切［4］，许多学者试图通过对免疫相关基因（immune-related genes，IRGs）或免疫相关lncRNA（immune-related lncRNAs，IRlncRNAs）的生物信息学分析来评估病人风险及预后，如Hong等［5］构建了12个IRlncRNAs的HCC预后模型，刘航睿等［6］构建了7个IRGs的食管癌预后模型，这些由IRGs或转录产物构成的模型在病人预后评价中具有良好的应用价值，但这些模型多只涉及一种类型的RNA，多类型RNA的联合预测模型研究较少。考虑到HCC免疫学研究的重要性，以及癌症发生发展中不同分子间的密切协同作用［7-8］，本研究首次构建了IRGs及IRlncRNAs联合的HCC预后模型，以期为HCC免疫学机制研究及预后预测提供一定参考。

1 资料与方法

1.1 数据获取及处理 2022年6—8月通过TCGA官网下载HCC转录组及与对应的临床数据；https：//www.immport.org/网站下载IRGs；根据Ensembl基因转移格式文件将基因表达矩阵区分为mRNA及LncRNA表达阵列，Perl软件提取mRNA阵列中IRGs表达阵列；剔除生存信息缺失及存活天数小于30的病人；R软件“WGCNA”函数对生存相关的核心免疫基因聚类，筛选与预后显著相关的基因模块（模块中基因数＞40）；上述得到的核心基因纳入单因素Cox回归筛选出预后相关的IRGs（P＜0.001为过滤条件），在预后相关的IRGs和lncRNAs之间进行共表达分析，得到IRlncRNAs（IRGs与lncRNA之间Pearson系数绝对值大于0.4，P＜0.001认为是IRlncRNAs）。

1.2 风险评分模型的建立及验证 套索算法（LASSO）是一种变量选择方法，通过构造一个惩罚系数得到比较精练的模型。使用R软件“glmnet”包使用R“caret”包按1∶1的比例将总样本随机分为训练集和验证集，在训练集中使用R“glmnet”包对上述得到的IRGs和lncRNAs进行Lasso回归分析，创建预后模型。风险评分公式=∑Ni=1（Ei×Ci）（Ei为每个元素的表达量，Ci为相对应的回归系数）。根据评分中位值将所有样本分为高风险组和低风险组，分析高风险组与低风险组HCC病人临床病理特征的差异（χ2检验）；再使用R软件绘制训练集组、测试集组、总样本组病人Kaplan-Meier曲线、时间依赖性受试者操作特征（ROC）曲线、风险因子图对模型进行评价；最后将HCC病人临床病理特征及模型风险评分纳入多因素Cox回归分析，判断模型预测病人预后的价值（P＜0.05为差异有统计学意义）。

1.3 不同风险组生存差异的机制分析 为了解不同风险组病人预后差异的分子机制，使用R语言limma包筛选高低风险病人差异表达的基因（FDR＜0.05，│Log2FC│＞1为过滤条件），并对差异基因进行GO（基因本体论）和KEGG（京都基因与基因组百科全书）分析，绘制富集结果的柱状图和气泡图（GO分析筛选条件为FDR＜0.05，基因计数＞5；KEGG通路分析筛选条件FDR＜0.05）。

2 结果

2.1 预后相关的IRGs及IRlncRNAs 剔除预后信息缺失及生存天数＜30 d的样本后共得到342个有效样本进行WGCNA分析，结果显示Power值为3时为最优无尺度网络，动态剪切法得到6个模块（图1A），黄色模块及蓝色模块中的基因为预后相关的核心基因（图1B），模块中共提取388个预后相关的核心基因，上述基因与生存数据合并后使用R“survival”包和单因素Cox回归进行IRGs筛选；共得到42个预后相关的IRGs，共表达分析得到42个预后相关IRlncRNAs（图2）。42个预后相关的IRGs均为HCC病人预后危险基因；在42个预后相关IRlncRNAs中，F11-AS1、AC015908.3、TMEM220-AS1、AP001065.3、AC115619.1是HCC病人预后的保护IRlncRNAs，其余为预后危险IRlncRNAs。

图1 通过WGCNA筛选肝细胞癌（HCC）预后相关的基因模块：1A为WGCNA聚类树；1B为与HCC病人生存相关的基因模块

图2 HCC预后相关的IRGs（A）和IRlncRNAs（B）森林图

2.2 LASSO回归构建模型 在训练集中对上述预后相关的IRGs及IRlncRNAs进行LASSO回归分析，选取误差最小的λ值确定最佳模型（图3），筛选出与预后显著相关的6个IRGs（NFYC、PSMD1、PSMD6、BCL10、TXLNA、S100A16）和7个IRlncRNAs（MSC-AS1、AC015908.3、AP001065.3、AL139384.1、AL031985.3、ZFPM2-AS1、AL117336.2），计算各个元素的回归系数，构建预后模型。

图3 基于IRGs和IRlncRNAs构建的HCC预后模型LASSO回归构建模型：3A为LASSO回归中的调谐平行度；3B为特征相对于log2(1ambda)的LASSO系数剖面

风险评分=0.062×NFYC+0.127×PSMD1+0.061×PSMD6+0.024×BCL10+0.1×TXLNA+0.009×S100A16+0.031×MSC-AS1-0.334×AC015908.3-0.007×AP001065.3+0.04×AL139384.1+0.053×AL031985.3+0.083×ZFPM2-AS1+0.444×AL117336.2。在模型中，NFYC、PSMD1、PSMD6、BCL10、TXLNA、S100A16、MSC-AS1、AL139384.1、AL031985.3、ZFPM2-AS1、AL117336.2为病人风险评分的危险变量，AC015908.3、AP001065.3为病人风险评分的保护变量。

2.3 模型的评价 根据评分中位值将342位病人分为高、低风险组。临床相关性分析显示：不同风险组病人的组织学分级Grade（P=0.001）、临床分期（P=0.005）、T分期（P=0.010）差异有统计学意义（表1）；模型在训练集（171例）、测试集（171例）及总样本集中分别进行验证，生存分析一致表明，高风险组较低风险组病人总生存期明显降低（均P＜0.05）（图4A～4C）；1年时间依赖性ROC曲线提示，训练集、测试集、总样本集曲线下面积95%CI分别0.85（0.79，0.91）、0.81（0.73，0.89）、0.83（0.78，0.88），提示该模型在预测HCC预后表现良好；风险热图显示低风险组病人生存状态更好（图4D）；多因素分析显示模型评分、临床分期是HCC病人独立的预后因子（均P＜0.05），见表2～4。

表1 模型风险评分与HCC病人临床病理特征的相关性/例

表2 训练集预后的单因素及多因素回归分析

表3 测试集预后的单因素及多因素回归分析

表4 总样本预后的单因素及多因素回归分析

图4 肝细胞癌（HCC）模型的评价及验证：4A为训练集组生存曲线；4B为测试集组生存曲线；4C为总样本组生存曲线；4D为风险得分；4E为生存状态；4F为模型变量的表达热图

2.4 高低风险病人差异基因的GO、KEGG分析 基于模型的不同风险组共存在1 241个差异表达基因。GO分析显示：差异基因主要富集的功能为有丝分裂、染色体分离、DNA复制起点的结合、微管蛋白结合、花生四烯酸单加氧酶活性、胶原结合等（图5A）；KEGG分析显示：差异基因主要富集的通路有PI3K-AKT、细胞周期、人乳头瘤病毒感染、黏着斑、ECM-受体相互作用、细胞衰老、蛋白质消化吸收等（图5B）。

图5 肝细胞癌不同风险组差异基因的基因本体论（GO）（5A）及京都基因与基因组百科全书（KEGG）（5B）分析棒状图及气泡图

3 讨论

尽管HCC的诊断和治疗方式取得较大进步，但具有不同分子特征的病人其治疗及预后存在明显差异，因此区分不同风险病人尤为重要。构建预后模型，能有效评估病人潜在的临床状态及结局，并筛选出与HCC发生发展及预后相关的分子标记，在个性化治疗和预后预测方面彰显优势。随着技术的进步，免疫治疗给人类抗癌带来了希望，同时在个性化的治疗时代，关键生物标志物及基因组的研究可引领肿瘤治疗的发展。鉴于此，部分研究者开始通过对IRGs或转录产物的研究来评估病人风险，但既往研究通常只关注单一类型的RNA的表达，如mRNA、lncRNA等，考虑到生物调控的复杂性以及不同RNA分子在癌症发展和进展中的密切联系，我们首次构建了HCC多类型RNA的预后模型，并分析了不同风险病人的临床特征，验证结果显示该模型预测HCC病人生存性能良好；多因素Cox回归提示该模型预测性能相对于传统临床特征是独立的；最后我们探讨了高低风险病人预后差异的可能分子机制，以期对HCC的临床决策及深入研究提供一定帮助。

在本模型纳入的6个IRGs和7个IRlncRNAs中，有部分已被报道与HCC的发展和预后有关。如Tan等［9］发现PSMD1的高表达与HCC的不良预后显著相关；Cheng等［10］发现BCL10在HCC组织中存在高频点突变，但具体致病机制仍不清楚；S100A16则可以参与HCC病人的体液免疫，并影响病人预后［11］。至于IRlncRNA方面，MSC-AS1被发现可预测HCC病人的预后、免疫细胞浸润水平及免疫治疗反应［12］；Kou等［13］发现MSC-AS1的下调可抑制HCC细胞增殖、迁移、侵袭和促进HCC细胞凋亡；ZFPM2-AS1可通过靶向调节miR-653/GOLM1轴［14］、miR-576-3p/HIF-1α轴［15］等途径促进HCC细胞增殖、侵袭，并与HCC不良预后相关。另外值得注意的是，AC015908.3的肿瘤干细胞相关功能新近在HCC中被鉴定出来［16］；AL031985.3则在HCC免疫、自噬等功能中发挥一定作用［17］，提示这些lncRNAs可能与HCC的发生发展密切相关。而关于NFYC、PSMD6、TXLNA、AP001065.3、AL139384.1、AL117336.2在HCC中的作用研究则较少，具体功能有待研究进一步揭示。最后GO分析显示高低风险病人风险差异基因主要参与细胞遗传物质的复制等事件，KEGG分析提示这些基因与某些肿瘤相关通路密切联系，如“PI3K-AKT”“细胞周期”等，表明HCC病人预后差异涉及多途径、多靶点，为HCC的多靶点联合治疗提供了理论支持。

HCC已进入免疫治疗时代，筛选免疫相关的生物标志物用于预测HCC病人预后非常重要和必要［18-19］。本研究是首个IRGs及IRlncRNAs联合的HCC预后模型的研究，该模型风险评分能够作为HCC病人的独立预后因素，将有助于临床医生对HCC病人的风险管理和决策制定。然而，我们的研究仍有一些不足。首先，本研究是基于国外数据库的分析，基于东西方HCC病人特征的不同，该模型在中国人群的应用是否同样具有价值尚未可知；其次，本模型是在公共数据的回顾性分析基础上构建，可能存在选择性偏倚；最后HCC具有高度异质性及病人个体差异，我们关注的队列是否具有代表性尚需在前瞻的独立临床队列中进行验证。总之，本研究表明，由IRGs及IRlncRNAs联合构建的预后模型在预测HCC病人的生存及临床特征上有一定价值，可能为HCC的免疫学研究及预后预测提供有用工具。