基于免疫相关IncRNA建立晚期头颈部鳞癌预后模型

吴嘉雯，唐加山

(南京邮电大学理学院，江苏 南京 210023)

头颈部癌症起源于上呼吸道和消化道的各种解 剖部位，居世界恶性肿瘤第六位。2020年全球范围 内头颈部癌症新增病例约93万人，死亡病例约47 万人[1]；其中，约90%的病例被归类为HNSCC。尽管近年来放疗、手术和辅助化疗等治疗策略取得了 进展，但是大约有三分之二的HNSCC患者在诊断 时已经处于局部晚期(III期和IV期)，预后较差[2]。 因此，建立可以准确预测晚期HNSCC患者预后的模型，对指导临床诊断和治疗有重要意义。

IncRNA被定义为大于等于200个核昔酸的RNA,目前还没有证据表明它们可以被翻译成肽。随着研究的深入,IncRNA的异常已经被证明具有抑制或促进肿瘤的作用，在肿瘤发展中发挥着不可或缺的作用[3-4]。此外,最近的研究表明，IncRNA在抗原呈递、免疫激活和免疫细胞浸润等癌症免疫的不同阶段具有重要作用[5]，因此免疫相关IncRNA引起了相当大的关注。免疫相关IncRNA作为新兴的癌症生物标志物已被用于多种癌症的诊断和生存预测,例如肝癌[6]、肺癌[7]、胃癌[8]等。同样，当前已有相关研究表明了免疫相关IncRNA在HNSCC中具有重要的预后价值[9-10]。然而，目前缺少关于晚期HNSCC患者预后相关的免疫IncRNA的研究。本文旨在识别与晚期HNSCC患者预后相关的免疫IncRNA,并基于这些IncRNA建立预后模型以改善晚期HNSCC的预后预测。

1 材料与方法

1.1 材料HNSCC患者的RNA测序(RNA Sequencing,RNA-seq)数据来自癌症基因组图谱(The Cancer Genome Adas,TCGA；https://portal,gdc.canc-er.gov/)中的HNSCC项目组，并从中提取了IncRNA表达数据;HNSCC患者的临床病理信息来自UCSC Xena(https://xenabrowser.net/)。将HNSCC患者的IncRNA表达数据和临床数据进行整理和合并，使得整理后的样本均有对应的表达数据和临床数据，并删除生存信息缺失和生存时间小于30天的样本。最后根据临床信息提取出临床III期和IV期的HNSCC样本。

1.2 方法

1.2.1 免疫相关IncRNA的提取 本文选取edgeR方法来对IncRNA进行差异表达分析,R软件中的edgeR是基于负二项分布的统计方法[11]，根据样本的IncRNA表达量，选择FDR<0.05和|log2FC|≥1作为阈值筛选出在肿瘤样本和非肿瘤样本之间有明显差异表达的IncRNA。另外从免疫学数据库(Immunology Database and Analysis Portal,ImmPort)(https：//www.immport.org/home)网站获得免疫相关基因列表，进行免疫基因与IncRNA的Pearson相关性分析，以相关系数|R|>0.4且P<0.001为筛选条件得到HNSCC免疫相关IncRNA。

1.2.2 关键预后免疫相关IncRNA的筛选 首先在训练集中对上述确定的免疫相关IncRNA进行单因素Cox回归分析，以P值小于0.05为标准确定与晚期HNSCC患者总生存期(overall survival,OS)相关的免疫IncRNA。然后基于这些预后相关的免疫IncRNA建立Lasso-Cox回归模型,并进行10折交叉验证确定最优模型,进一步筛选出与晚期HNSCC患者预后密切相关的免疫IncRNA。同样，在训练集中对确定的免疫相关IncRNA进行Coxboost分析来选择与晚期HNSCC患者预后相关的免疫IncRNA。基于R语言中的“Coxboost”包构建Coxboost模型，并使用cv.CoxBoost函数进行5折交叉验证来选择最优提升步数。最后比较筛选得到的两组与晚期HNSCC患者预后相关的免疫IncRNA,确定用于建立预后模型的关键免疫IncRNA。

1.2.3 预后模型的构建 基于筛选出来的关键预后免疫相关IncRNA,利用随机生存森林算法建立晚期HNSCC患者的预后模型,并计算每个患者的预后风险值。绘制模型的5年OS的时间依赖性受试者工作特征曲线(time depesndent receiver operating characteristic curve,timeROC)曲线，计算使Youden指数最大时的风险值作为阈值，并根据风险值的阈值将训练集和测试集中的患者分为高风险组(风险值大于阈值)和低风险组(风险值小于等于阈值)。

1.3 统计学分析本文所有的统计分析和可视化均基于R软件(4.1.2版)进行,P值小于0.05被认为具有统计学意义,P值小于0.01认为具有显著差异。

2 结果

2.1 数据处理从UCSC Xena网站下载了612名HNSCC患者的临床信息，并将临床信息和RNA-seq数据进行匹配，得到539个匹配样本;其中包括495个肿瘤样本和44个正常组织样本。接着去除生存信息缺失和生存时间小于30 d的样本，并从中提取III期和IV期样本，最终得到362个晚期HNSCC样本。

根据基因注释文件，从HNSCC的RNA-seq数据中得到15 878个IncRNA的表达数据，通过差异表达分析(FDR<0.05且|log2FC|≥1)得到1 729个差异表达的IncRNA。根据ImmPort数据库下载的免疫相关基因信息，从基因表达数据中提取了1 279个免疫相关基因的表达数据。采用Pearson相关分析筛选与免疫相关基因相关的IncRNA,结果确定了988个免疫相关IncRNA(|R|>0.4,P<0.001)。

2.2 关键预后相关免疫IncRNA的筛选及预后模型的构建按照7∶3的比例将样本随机分为训练集和测试集。在训练集中，我们对包含988个免疫相关IncRNA的生存数据进行了单因素Cox回归分析，以确定与预后相关的免疫IncRNA；根据Wald检验的P值小于0.05为标准，筛选与预后相关的免疫IncRNA,结果得到159个免疫相关IncRNA。再基于这159个免疫相关IncRNA建立Lasso-Cox回归模型，通过10折交叉验证的Lasso-Cox回归分析结果如图1所示。图1A中的横坐标为A的对数值，纵坐标为方差值，最上方为Lasso-Cox回归分析后的IncRNA的数量，选择使模型达到最小损失时的惩罚参数λ(λ=0.118)，得到8个系数非零的预后相关免疫IncRNA,分别是LINC01305、RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-30K9.5、RP11-116D17.3和RP11-445F12.1。

图1 Lasso-Cox回归分析筛选关键预后相关免疫IncRNA

同样，在训练集上建立Coxboost模型，进行10次提升之后，得到7个系数非零的预后相关免疫lncRNA,为LINC01305、RP11-30P6.6、RP11-65M17.3、RP11-497E19.2、RP11-30K9.5、RP11-116D17.3和RP11-445F12.1。

对筛选得到的两组预后相关免疫IncRNA进行比较,可以看到基于Cox回归和Lasso-Cox回归筛选出来的与预后相关的免疫IncRNA集合中仅比基于Coxboost分析筛选得到的预后相关的免疫lncRNA集合多了IncRNA RP11-890B15.2。表明了UNC01305、RP11-30P6.6、RP11-65M17.3、RP11- 497E19.2、RP11-30K9.5、RP11-116D17.3和RP11-445F12.1这7个免疫相关IncRNA与预后具有密切的联系。为了确定RP11-890B15.2对于预后的重要性，我们基于RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-116D17.3、RP11-445F12.1、LINC01305和RP11-30K9.5这8个IncRNA,利用随机生存森林构建预后模型，并根据VIMP值对于这8个预后相关免疫IncRNA进行重要性排序。通过网格搜索计算并选择选择使随机生存森林的袋外错误率最小的mtry和node size组合作为该模型的最优参数。根据计算结果可以得到，在mtry=45,node size=1时构建的随机生存森林模型的袋外错误率达到最低，此时的袋外错误率为33.3%,并且模型的袋外错误率在1 000棵生存树时已经趋于稳定(见图2A)。计算变量VIMP值并根据VIMP值对IncRNA的重要性进行排序，从图2B和表1中可以看出这8个免疫相关IncRNA对于预后模型的重要性均是正值，提示了这8个免疫相关IncRNA对于患者预后都具有一定的影响，因此将这8个免疫相关IncRNA全部纳入预后模型。

图2 随机生存森林构建预后模型

表1 8个免疫相关IncRNA的重要性值

2.3 预后模型的验证在训练集上,根据上述基于随机生存森林建立的预后模型计算每个患者的风险值，并根据风险值的阈值(cutoff=33.07)将训练集中风险值大于阈值的患者归为高风险组，反之则归为低风险组;训练集和测试集中患者风险值的分布、相应的生存状态和8个免疫相关IncRNA的表达水平的热图分别如图3和图4所示，从图中可以看出高风险组发生死亡的人数多于低风险组，且高风险组患者的生存时间少于低风险组。热图显示,UNC01305和RP11-30K9.5在低风险组中高表达,而RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-116D17.3和RP11-445F12.1在高风险组中高表达。Kaplan-Meier生存分析结果显示，低风险组患者总生存率显著高于高风险组患者(P<0.001)(图5A)。在测试集中进行了相同的分析，并且得到了类似的结果(结果见图5B)。另外,在总数据集中，低风险组患者的中位生存期为2 900 d而高风险组患者的中位生存期882 d，可见高风险组相比低风险组的总体生存较差(见图5C)。

图3 训练集风险值分析

图4 测试集风险值分析

图5 高风险组和低风险组患者的Kaplan Meier生存曲线

另外，我们绘制了timeROC曲线并根据AUC值来评估模型的预后价值，计算得到训练集上1年、2年和3年OS的ROC曲线的AUC值分别为0.75、0.76和0.77(图6A)；测试集上1年、2年和3年ROC曲线的AUC值分别为0.61、0.66和0.64(图6B)；总数据集上1年、2年和3年ROC曲线的AUC值分别为0.71、0.73和0.73(图6C)。

图6 timeROC曲线

t-分布随机邻域嵌入(t-SNE)算法是一种降维技术,经常用于将高维数据映射到二维空间从而使其可视化。我们使用t-SNE探索了训练集(图7A)和测试集(图7B)中高风险组和低风险组患者的分布,从图中可以直观地看到，具有不同风险的患者被很好地分为两个集群。

图7 训练集和测试集的t分布随机邻域嵌入分析

2.4 预后模型在不同亚组中的验证为了进一步评估所构建的预后模型的可靠性，我们在整个数据 集中进行了分层分析。根据年龄小于70岁和大于等于70岁划分为不同的年龄组;分别在不同的年龄组内对局低风险组患者进行Kaplan-Meier生存分析并利用Log-rank检验进行高低风险组之间生存率的比较。根据Kaplan-Meier曲线可知,在小于70岁的分组中高风险组患者生存率显著低于低风险组患者的生存率(P<0.0001,结果见图8A)；根据Kap-lan-Meier方法估算得到的高风险组患者的中位生存期为2.5年,低风险组患者中位生存期为7.9年;可知低风险组患者的预后较高风险组患者更好。在年龄大于等于70岁的分组中,高风险组患者的中位生存期为2.3年，低风险组中位生存期为5.7年，高风险组患者生存率低于低风险组患者生存率，且具有统计学意义(P=0.011,见图8B)。同样,在不同性别组内进行Kaplan-Meier生存分析，根据分析结果可知，男性患者组(P<0.000 1,见图8C)和女性患者组(P=0.034,见图8D)的高、低风险组之间的生存率存在差异且具有统计学意义;其中在男性患者组中,高风险组患者的中位中位生存期为2.6年，低风险组患者中位生存期为7.9年;女性患者组中的高风险组患者的中位生存期为2.3年。综上可知,该预后模型在不同的年龄亚组和性别亚组中均具有良好的预测能力，表明了该模型是可靠的。

3 讨论

当前的数据表明，近年来HNSCC的发病率和死亡率在不断地上升，尽管诊断和治疗手段的进步使得早期患者得到了改善，但是大部分患者在确诊时处于晚期[12]，耐药性难以克服，使其预后不良。近年来,lncRNA成为肿瘤研究的热点，人们发现lncRNA在肿瘤的发生和发展中起着重要作用，lncRNA的异常表达被认为是影响肿瘤活性的关键因素;其中与免疫相关的IncRNA被证明与多种恶性肿瘤的预后有关[13-14]。然而，目前缺少关于晚期HNSCC患者预后相关的免疫IncRNA的研究。因此，本文的研究目的是通过对TCGA数据库中HNSCC晚期患者的相关数据进行分析，确定与晚期HNSCC预后相关的免疫IncRNA,用于晚期HNSCC患者的预后风险评估。

本研究中，使用单因素Cox回归、Lasso-Cox回归和Coxboost分析,确定了8个与晚期HNSCC预后相关的免疫IncRNA(LINC01305、RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-30K9.5、RP11-116D17.3和RP11-445F12.1),基于8个预后相关免疫IncRNA利用随机生存森林建立了晚期HNSCC患者的生存预后风险模型,并根据风险值的阈值将晚期HNSCC患者划分为高风险组和低风险组。根据8个IncRNA的表达量热图可知，IncRNA LINC01305和RP11-30K9.5在低风险组中高表达,提示IncRNA LINC01305和RP11-30K9.5的高表达量可能与患者较好的预后有关;而RP11-30P6.6、RP11-890B15.2、RP11-65M17.3、RP11-497E19.2、RP11-116D17.3和RP11-445F12.1在高风险组中高表达，提示它们的高表达量可能与患者较差的预后有关。根据Kaplan-Meier曲线可以看到，高风险组患者的生存率显著低于低风险组患者的生存率(P<0.01),其中根据基于整个数据集的Kaplan-Meier分析结果可知，高风险组患者的3年生存率约为41%,中位生存期约为2.4年；而低风险组患者的3年生存率约为85%,中位生存期约为7.9年。高风险组患者的3年生存率比低风险组患者的3年生存率高了44%,可见高风险组的患者预后较低风险组更差。同时，在训练集、测试集和总数据集中，模型3年OS的timeROC曲线的AUC值分别为0.77、0.64和0.73,表明该模型具有较好的预测能力。为了进一步验证预后风险模型在不同临床特征分层中的适用性，分别在不同性别组和不同年龄组中对高低风险组患者的生存率进行了分析，分析结果表明,基于8个免疫相关IncRNA的预后风险模型不受性别和年龄因素的影响，在不同的分层中的高低风险组患者生存率的差异均具有显著性。

另外通过查阅文献发现,IncRNA LINC01305可以促进宫颈癌的进展，且LINC01305的高表达与宫颈癌患者的低生存率相关[15]；且有研究证明了UNC01305可以促进食管鳞状细胞癌的转移和增殖。迄今为止还没有相关报告证明这8个免疫相关IncRNA与HNSCC的关系，它们在HNSCC中的作用机制还需要进行进一步的探索。

这项研究也存在一定的局限性。首先,本文仅基于TCGA-HNSCC数据进行分析，样本量较少可能存在分析偏差;另外,这是一项纯粹的统计分析和数据挖掘研究，构建的预后风险评估模型在临床预测中是否可行还有待进一步验证。

综上所述，本文筛选出8个与晚期HNSCC患者预后相关的免疫IncRNA，并基于这8个免疫相关IncRNA成功构建了预后模型,经验证表明该预后模型能较好的预测晚期HNSCC患者的预后情况。可能为未来的研究和临床实践提供新见解。