mRNAs列线图模型对宫颈腺癌患者生存状况的预测价值

张钗红，关德凤，杨永秀

宫颈癌是女性最常见的癌症之一，发病率和病死率位居女性恶性肿瘤第4位[1]，仍是成年女性生命安全的重大威胁。近年来，随着宫颈癌规范化筛查的开展和疫苗的应用，宫颈鳞癌的发病率不断下降，但宫颈腺癌(cervical adenocarcinoma，CAC)发病率仍不断上升，占宫颈癌的10%～20%[2-3]，并且其与人乳头瘤病毒的关系尚有争议[4-5]，发病人群呈年轻化趋势，早期筛查困难，具有较鳞癌更高的卵巢转移及远处转移率，宫颈腺癌组织学已成为宫颈癌预后不良的因素[6-7]，针对腺癌组织学特异性的基因靶向疗法是宫颈癌新的研究热点，也是改善患者生存期新的切入点。本研究构建列线图预后模型，使用GSEA_4.0.1软件进行模型mRNAs在CAC中的癌症和效应特征基因集富集分析(gene set enrichment analysis, GSEA)，探索模型mRNAs在CAC中的作用机制及模型mRNAs之间相互关系，寻找稳定可靠的CAC靶向治疗和预后评估的关键因子，报道如下。

1 资料与方法

1.1 数据准备 本研究于2019年8—10月在兰州大学第一医院实施，使用R3.6.1软件分析癌症基因组图谱数据库(the cancer genome atlas，TCGA)中宫颈腺癌mRNA表达数据。在TCGA官网(https://www.cancer.gov/tcga)获取CAC的mRNA表达数据建立标准化表达矩阵并加载R3.6.1软件edgeR包，采用双侧t检验，筛选显著差异表达mRNAs；获取CAC患者临床资料，合并mRNAs差异表达资料构建预后预测模型。

1.2 筛选预后相关mRNAs 运行R3.6.1软件survival包、glmnet包进行单因素Cox回归、Lasso回归分析、多因素Cox逐步回归分析，进一步筛选预后相关mRNAs，使用survminer包进行Log-rank检验筛选CAC患者独立的预后影响因子，计算C-指数评价模型。

1.3 列线图模型构建及评价 在R3.6.1软件中加载rms包和survival包构建CAC患者3年和5年生存率的列线图预测模型。加载timeROC包和survival包依据46个CAC样本的风险评分，计算3年和5年生存率的AUC值。

1.4 模型mRNAs分析 在R3.6.1软件中加载limma包和beeswarm包采用Wilcox检验对模型mRNAs在正常组和肿瘤组中的标准化表达数据进行分析。加载survminer包和survival包，使用中位数法将46个CAC样本的风险评分分为高风险组和低风险组，进行高、低风险组生存分析。加载survival包，使用中位数法分别将5个mRNAs分为高、低表达组，进行mRNAs高、低表达组生存分析。构建模型mRNAs的GSEA输入文件：表型文件和表达谱文件，使用GSEA_4.0.1软件设置基因集数据为c6.all.v7.0.symbols.gmt和h.all.v7.0.symbols.gmt进行高、低表达组的癌症、效应、免疫特征基因集富集分析，使用R3.6.1软件plyr、ggplot2和grid包绘制多GSEA富集图。

2 结 果

2.1 mRNAs表达数据和临床资料 在TCGA数据库下载CAC的mRNAs表达谱资料，包含3个正常样本和48个CAC样本，加载R3.6.1软件进行差异表达分析得到CAC中1 856个|logFC|>2且FDR<0.01的显著差异mRNAs，即649个显著上调mRNAs和1 207个显著下调mRNAs。下载CAC患者临床资料，剔除生存时间为0的2组样本，共入组46个CAC样本，用于模型构建。

2.2 筛选CAC预后相关mRNAs 使用R软件对CAC的mRNAs表达与生存数据矩阵进行单因素Cox回归分析，筛选得到P<0.01的22个CAC预后相关mRNAs，提取22个mRNAs的表达生存数据进行Lasso回归分析，得到最小Lambda为0.110 273 7，包含10个CAC预后相关mRNAs，分别是ADAMTSL4、ANKRD53、CDC25A、CILP、CPQ、CTSF、DST、GLIPR2、MYH13、RNF150。提取这10个mRNAs的表达生存数据进行多因素Cox逐步回归分析，获得5个CAC预后相关mRNAs及其在46个CAC样本中的风险评分，该模型风险评分=0.944 74×CILP+3.866 33×MYH13-3.776 3×ANKRD53-1.912 99×GLIPR2-0.641 92×CPQ，经Log-rank检验P<0.01，其中ANKRD53、CILP、GLIPR2、MYH13可作为CAC患者预后评估的独立影响因素，P值分别为0.009、0.024、0.007和0.005，而CPQ的P值为0.131，不能作为患者独立预后因子。

2.3 CAC患者生存预测模型及其评价 使用R软件构建CAC患者3年生存率和5年生存率的列线图预测模型(图1)。依据46个CAC样本的风险评分计算模型C-指数为0.95，3年生存率的AUC为0.950 1，5年生存率的AUC为0.949 7，该模型预测效果较好。

2.4 模型mRNAs分析

2.4.1 模型mRNAs的TCGA表达分析:从CAC的mRNA标准化表达数据矩阵中提取模型mRNAs的表达数据，使用R3.6.1软件进行正常组和肿瘤组的表达差异分析显示，2组中ANKRD53、CILP、GLIPR2、CPQ、MYH13的表达差异均有统计学意义(P均<0.01)。

2.4.2 风险评分与生存预后的关系：使用R软件进行高、低风险组生存分析显示，低风险组患者生存预后显著优于高风险组(P<0.01)。

2.4.3 模型mRNAs表达情况与患者预后的关系：使用R软件进行mRNAs高、低表达组生存分析显示，ANKRD53、CLIP、GLIPR2、CPQ高表达组和MYH13低表达组提示CAC患者预后良好，其中ANKRD53、CLIP、GLIPR2高、低表达组差异有统计学意义(P=0.023、0.028、0.030)。

2.4.4 模型mRNAs的基因集富集分析

2.4.4.1 癌症特征基因集分析(oncogenic signatures gene sets) 使用GSEA_4.0.1软件输入模型mRNAs的表型文件和表达谱进行癌症基因集分析显示，5个模型mRNAs共存在33个共同显著相关(P<0.05)的癌症特征基因集，其中VEGF_A_UP.V1_UP基因集与4个模型mRNAs均显著相关(图2)，RB_P107_DN.V1_DN、JAK2_DN.V1_UP、CYCLIN_D1_KE_.V1_DN、JNK_DN.V1_DN、P53_DN.V1_DN、BMI1_DN_MEL18_DN.V1_DN、KRAS.KIDNEY_UP.V1_UP、PTEN_DN.V1_UP基因集分别与3个模型mRNAs显著相关，24个基因集分别与2个模型mRNAs显著相关。

2.4.4.2 效应特征基因集分析(hallmark gene sets) 使用GSEA_4.0.1软件对5个mRNAs的表型文件和表达谱进行效应特征基因集分析显示，5个模型mRNAs共存在9个共同显著相关(P<0.05)的效应特征基因集，其中MYC TARGETS V1和MTORC1 SIGNALING分别与3个模型mRNAs显著相关，MYC TARGETS V2、ANGIOGENESIS、EPITHELIAL MESENCHYMAL TRANSITION、HEDGEHOG SIGNALING、UNFOLDED PROTEIN RESPONSE、FATTY ACID METABOLISM、MYOGENESIS分别与2个模型mRNAs显著相关(图3)。

3 讨 论

近年来，宫颈腺癌发病率不断上升，预后不良，复发率、转移率、病死率均较高[6-9]，有效的预后风险评估、精准的基因靶向治疗方案、提高治愈率、延长生存期仍是现阶段宫颈腺癌研究的重点和热点。

现阶段，宫颈腺癌患者的预后评估工作主要基于不同的临床病理参数，如是否绝经、子宫旁受累、脉管内瘤栓、淋巴结转移、Silva分型系统等[6, 10]，不同的影像学参数，如磁共振弥散加权成像和波谱分析等[11-12]。本研究通过整合TCGA中46例CAC患者的mRNAs表达和生存数据，使用R软件进行单因素Cox、Lasso、多因素Cox逐步回归分析，构建预后评估准确度较高的分子模型，经模型评价预测准确性高，C-指数为0.95，3年生存率的AUC为0.950 1，5年生存率的AUC为0.949 7，其中ANKRD53、CILP、GLIPR2、MYH13可作为CAC患者预后评估的独立影响因素，可初步筛选预后不良需进一步行预防性靶向干预的患者。

目前，宫颈腺癌的分子靶向治疗研究主要集中在探索可以作为患者预后标志物和治疗靶标的潜在基因。近年来，高通量测序技术已广泛应用于CAC的基因分析，TCGA数据库31例CAC患者的体细胞突变数据显示,90.32%的CAC组织样本存在体细胞突变，涉及17 342个基因的8 867个染色体突变。基因突变与CAC患者预后密切相关。Jiang等[13]发现，宫颈腺癌较鳞癌具有更高的KRAS非同义突变率和更短的无复发生存期。精准靶向性干预异常基因可有效改善CAC治疗效果。Takiguchi等[14]进行蛋白质印迹分析发现，Src基因在CAC细胞系(Hela和TCO-2)中高表达，使用比色法测定活细胞计数发现，Src抑制剂联合化疗药可显著抑制CAC细胞生长增殖。Akimoto等[15]也发现claudin-1在宫颈腺癌呈高表达，敲除claudin-1可诱导细胞凋亡，显著抑制CAC细胞的增殖、迁移和侵袭。本研究通过分析高通量测序得到的CAC全基因组表达网络，筛选出高度异常表达的mRNAs谱，结合患者生存数据，获取与CAC患者预后高度相关的新型分子标记构建生存预测模型，并通过对比模型mRNAs与全基因组的表达数据，进行模型mRNAs的癌症和效应特征基因富集分析，探索模型mRNAs在CAC中的作用机制及模型mRNAs之间的相互关系。癌症特征基因集分析发现，5个模型mRNAs共存在33个共同显著相关的癌症特征基因集，其中VEGF_A_UP.V1_UP基因集与4个模型mRNAs均显著相关，8个基因集分别与3个模型mRNAs显著相关，24个基因集分别与2个模型mRNAs显著相关。ANKRD53、CILP、CPQ和GLIPR2等4个模型mRNAs在癌症特征方面的联系更加紧密，与VEGF_A_UP.V1_UP基因集均显著相关。MYH13主要通过CPQ与ANKRD53、CILP和GLIPR2发生联系。效应特征基因集分析发现，5个模型mRNAs共存在9个共同显著相关的效应特征基因集，其中MYC TARGETS V1和MTORC1 SIGNALING分别与3个模型mRNAs显著相关，7个基因集分别与2个模型mRNAs显著相关。ANKRD53、CILP、CPQ和MYH13之间联系紧密，共存在8个相关基因集。GLIPR2主要通过CPQ与ANKRD53、CILP和MYH13发生联系。总之，癌症和效应特征基因集富集分析显示，5个模型mRNAs联系紧密，主要可能通过VEGF_A_UP.V1_UP、MYC TARGETS V1和MTORC1 SIGNALING等途径在CAC中发挥作用并影响患者预后。

综上所述，本研究通过分析TCGA中CAC患者的mRNAs表达和临床数据集，构建预后评估准确度较高的列线图模型，用于筛选预后不良需进一步干预的患者，结合癌症和效应特征基因集富集分析，探索模型mRNAs在CAC中的作用机制及模型mRNA之间的相互关系，从而为CAC患者精准靶向性个体化治疗和预后评估提供新的线索。

(致谢：大连医科大学基础医学院全昱冲硕士在实施R软件分析中的指导，谨此致谢)

利益冲突：所有作者声明无利益冲突

作者贡献声明

张钗红：课题设计、实施，数据分析，论文撰写；关德凤：资料收集、数据核对；杨永秀：论文审核