基于公共数据库构建甲状腺癌差异基因Cox比例风险回归模型对患者预后评估价值

姜雪丽，赵 辛，赵铁铮

甲状腺癌是一种好发于中老年女性的内分泌系统常见恶性肿瘤，多数病理类型为甲状腺乳头状癌[1-4]。据统计，2017年全球新发甲状腺癌病例占所有恶性肿瘤新发病例3.4%，在女性恶性肿瘤新发病例中甲状腺癌居第5位[5]。我国每年甲状腺癌新发病例和死亡病例数占全世界15%左右[6]。近年甲状腺癌患者病死率逐年增加[7-9]。目前，临床上对甲状腺癌患者预后判断主要基于临床病理特征，但由于恶性肿瘤有明显异质性，且患者个体差异较大，故此预测手段具有明显局限性[10]。近年，肿瘤基因和蛋白表达预测患者预后正逐步受到重视。

本研究通过对GEO和TCGA数据库中甲状腺癌患者转录组数据进行联合分析，筛选出甲状腺癌中差异表达的基因;随后通过单因素Cox回归分析，筛选出与预后相关的差异基因;最后通过多因素Cox回归分析成功构建5个差异基因预测甲状腺癌患者预后的模型，并结合患者临床资料对模型进行验证，用于判断甲状腺癌患者预后，现报告如下。

1 资料与方法

1.1资料收集和差异分析 在GEO数据库中以“Thyroid cancer”为关键词，选取人甲状腺癌组织样本转录组数据集GSE138198和GSE50901，使用R软件及Sva功能包校正批次，利用Limma包分析差异基因，认为符合|log2(Fold Change)|>0.585且P<0.05的基因为差异基因。在TCGA数据库中下载人甲状腺癌RNA-Seq数据，共含有58份正常甲状腺样本和509份甲状腺癌样本，利用R软件及Limma包分析差异基因，认为符合|log2(Fold Change)|>2且P<0.01的基因为差异基因。甲状腺癌患者的临床资料利用R软件进行整理，包括性别、年龄、临床分期、风险值评分、生存时间和生存状态等，去除临床资料不完整患者，最终保留498例采用单因素和多因素Cox回归分析进行甲状腺癌患者临床病理特征与预后关系的研究。

1.2Cox比例风险回归模型建立和验证 利用R软件将TCGA数据库中甲状腺癌患者随机分为两组，即试验组和验证组，使用Survival包对在GEO和TCGA数据库中均发生改变基因与患者预后的关系进行单因素Cox回归分析；使用Glmnet包对单因素Cox回归分析结果有差异的基因(P<0.05)进行LASSO回归分析；使用Survival和Survminer包对经LASSO回归分析筛选后的基因进行多因素Cox回归分析，进一步检验具有风险预测能力基因并计算各基因回归系数，从而构建判断患者预后的Cox比例风险回归模型。风险值为样本中各基因表达量与回归系数乘积之和，计算试验组、验证组及整体组(TCGA数据整体)各样本风险值，选取试验组风险值中位值为临界值，将各组分别分为高风险组和低风险组，绘制Kaplan-meier生存曲线，采用受试者工作特征(ROC)曲线评估模型预测甲状腺癌患者生存率准确性，借助R软件及Ggplot2、Survminer、SurvivalROC、Pheatmap等功能包进行数据可视化。

1.3GO和KEGG功能富集分析 将TCGA数据库中甲状腺癌基因表达矩阵按照风险值高低分为高风险组和低风险组两组，使用R软件及Limma包行差异基因筛选，符合|log2(Fold Change)|>0.585且P<0.05的基因为差异基因，利用R软件及Enrichplot、Org.Hs.eg.db、Cluster Profiler、Ggplot2等功能包进行GO和KEGG功能富集分析及数据可视化。

1.4统计学方法 应用R 3.6.1软件及相关功能包对所有数据进行处理分析，α=0.05为检验水准。

2 结果

2.1差异基因筛选 对GEO数据库中甲状腺癌数据集GSE138198和GSE50901进行批次校正后发现，与正常甲状腺组织相比，甲状腺癌组织中有897个基因表达发生改变，其中506个基因上调，391个基因下调，见图1a；结合TCGA数据库分析发现，在此两个数据库甲状腺癌组织中有241个基因均上调，207个基因均下调，见图1b。

图1 甲状腺癌组织和正常甲状腺组织差异基因筛选

2.2差异基因单因素Cox回归分析 对试验组差异基因进行单因素Cox回归分析显示，57个差异基因与甲状腺癌患者预后相关(P<0.05)，其中27个差异基因为高风险基因，30个差异基因为低风险基因，见图2。

图2 甲状腺癌患者差异基因单因素Cox回归分析

2.3Cox比例风险回归模型构建 试验组单因素Cox回归分析中与甲状腺癌患者预后有关的差异基因经LASSO回归分析筛选出6个候选差异基因，对其进行多因素Cox比例风险回归模型构建，最后发现PHLDA2、GPR137B、PORCN、MAPK4和TSPYL2共5个基因参与模型构建，其中PHLDA2为低风险基因，GPR137B、PORCN、MAPK4和TSPYL2为高风险基因,见表1。风险值=(PHLDA2×-0.1028)+(GPR137B×0.0880)+(PORCN×0.1112)+(MAPK4×0.2403)+(TSPYL2×0.1465)。基于风险值将试验组分为高风险组和低风险组两组，生存分析发现，高风险组生存时间和生存率低于低风险组，见图3a，且此模型中试验组1年生存率的ROC曲线下面积为0.929，见图3b。

表1 甲状腺癌患者差异基因Cox比例风险回归模型

图3 甲状腺癌患者差异基因Cox比例风险回归模型构建分析3a和3b均为试验组

2.4Cox比例风险回归模型验证 利用验证组及整体组对Cox比例风险回归模型进行验证，生存分析发现，验证组和整体组高风险组生存时间和生存率低于低风险组，见图4a和4b；此外，验证组和整体组1年生存率的ROC曲线下面积分别为0.680和0.773，见图4c和4d。

图4 甲状腺癌患者差异基因Cox比例风险回归模型验证4a和4c为验证组，4b和4d为整体组

2.5GO和KEGG功能富集分析 对TCGA数据库高和低风险组进行差异基因分析发现，152个基因在高和低风险组甲状腺癌样本中差异表达。差异基因GO功能富集分析发现，受体介导的内吞作用和先天免疫应答激活信号传导等显著富集，见图5a；KEGG功能富集分析发现，mTOR信号通路、甲状腺激素信号通路和HIF1信号通路等显著富集，见图5b。

图5 甲状腺癌样本差异基因GO和KEGG功能富集分析5a为GO功能富集分析，5b为KEGG功能富集分析

2.6甲状腺癌患者临床病理特征与预后关系 采用单因素和多因素Cox回归分析评估甲状腺癌患者临床病理特征与预后关系。单因素Cox回归分析结果显示，甲状腺癌患者年龄、临床分期、风险值评分与预后有关(P<0.01)；多因素Cox回归分析结果显示，甲状腺癌患者年龄和风险值评分与预后相关(P<0.05或P<0.01)，见表2。

表2 甲状腺癌患者临床病理特征与预后关系

3 讨论

目前，临床上对恶性肿瘤患者预后判断主要基于临床病理特征[11]，具有一定主观性，且患者预后个体差异较大，故构建科学预后评估模型对恶性肿瘤治疗效果评估具有重要意义[12]。Cox比例风险回归模型是以生存时间和生存状态为应变量，可同时分析多种因素对生存情况影响的一种半参数回归模型。近年来，越来越多的学者将其应用于恶性肿瘤患者预后预测研究中[13]，通过对肿瘤样本进行测序分析，结合患者临床病理特征，构建基于microRNA、LncRNA、mRNA和蛋白等的预后预测模型[14]，在临床诊疗中具有广阔应用前景，吸引了广大科研工作者的注意。

本研究通过对公共数据库中甲状腺癌样本的转录组数据进行联合分析，筛选出在甲状腺癌组织和正常甲状腺组织中差异表达的基因，对差异基因进行单因素Cox回归分析显示，有57个差异基因与甲状腺癌患者预后相关，其中27个差异基因为高风险基因，30个差异基因为低风险基因，随后通过LASSO回归分析和多因素Cox回归分析构建由PHLDA2、GPR137B、PORCN、MAPK4和TSPYL2共5个基因组成的Cox比例风险回归模型。基于此模型，将试验组、验证组和整体组分别分为高风险组和低风险组，生存分析发现各高风险组生存时间和生存率低于各低风险组，采用ROC曲线评估此模型预测甲状腺癌患者生存率的准确性，发现此模型具有较高准确性。多因素Cox回归分析结果显示，甲状腺癌患者年龄和风险值评分与预后相关。上述结果表明此模型具有一定临床应用价值，可用于甲状腺癌患者预后判断。此外，通过对高和低风险甲状腺癌样本进行差异基因筛选和功能富集分析发现，受体介导的内吞作用和先天免疫应答激活信号传导等显著富集，mTOR信号通路、甲状腺激素信号通路和HIF1信号通路等显著富集。既往研究也表明，mTOR信号通路和HIF1信号通路等与甲状腺癌细胞的增殖、迁移及侵袭等恶性行为密切相关[15]。

对本研究构建模型的5个差异基因进行文献检索发现，尽管其在甲状腺癌中的作用报道较少，但在其他类型肿瘤中有较多研究。GPR137B在本研究模型中是一个高风险基因，在甲状腺癌组织中高表达(数据未展示)，其在肿瘤中的作用以往未见报道，但敲低其同一家族基因GPR137已被证实可以抑制卵巢癌、胰腺癌和肝癌等细胞的增殖[16-17]；PORCN蛋白介导WNT蛋白棕榈酰化，对于WNT蛋白的分泌及WNT信号通路的激活具有重要意义，而WNT信号通路激活与甲状腺癌恶性进展密切相关[18]；MAPK4过表达与肺腺癌、膀胱癌、低级别胶质瘤和甲状腺癌生存差异相关，过表达MAPK4通过激活AKT/mTOR信号通路促进肿瘤恶性进展[19]。本研究结果与上述研究结果基本相符，提示GPR137B、PORCN和MAPK4在甲状腺癌中具有重要生理功能。然而，由于甲状腺癌研究队列的不足，本研究将TCGA数据库中甲状腺癌患者通过随机分组的方式构建模型并验证，故需要更多其他甲状腺癌队列研究来进一步验证本研究所构建的模型，且GPR137B、PORCN和MAPK4等基因对甲状腺癌细胞增殖、迁移和侵袭等的影响需要后续细胞和动物实验进一步探讨。

总之，本研究基于公共数据库成功构建5个差异基因的甲状腺癌Cox比例风险回归模型，具有较高的准确性和可靠性，有助于临床医生判断甲状腺癌患者预后。