基于生物信息学筛选与宫颈癌免疫相关的分子标志物

罗永金，胡晓霞，王 丹，刘 珍，张 侃，洪翠萍

宫颈癌是常见的妇科恶性肿瘤之一[1]，根据2018年的全球数据，每年约有31.1万患者死于宫颈癌[2-3]。高危型人乳头瘤病毒(human papilloma virus，HPV)的持续感染是宫颈癌发生的主要原因[4-6]。但除了HPV感染外，患者的遗传免疫因素也起到了重要作用。基因多态性的异常和免疫功能失调等因素与持续的HPV感染在宫颈癌的发生过程中具有协同作用[7-8]。肿瘤微环境(tumor microenvironment,TME)是肿瘤细胞所处的细胞环境，它由免疫细胞、间充质细胞、内皮细胞、炎性介质和细胞外基质(extracellular matrix,ECM)组成[9-10]。TME中的细胞和分子处于一个动态的过程，反映了癌症的进化本质，并共同促进肿瘤的免疫逃逸、生长和转移[11-12]。免疫细胞和基质细胞是两种主要类型的非肿瘤成分，被认为在肿瘤的诊断和预后中具有重要的价值[13]。随着生物学大数据时代的发展，通过生物学、计算机科学和信息技术的结合，对包含大量数据的生物信息进行分析成为可能。2013年，Yoshihara等[14]设计了一种ESTIMATE算法,利用表达数据估算恶性肿瘤组织中的基质细胞和免疫细胞。在此算法中，研究者通过计算免疫细胞和基质细胞中特异性分子生物标志物的表达特征，通过免疫和基质评分来预测TME。近年来，有报道将ESTIMATE算法应用于胶质母细胞瘤[15]、透明细胞肾细胞癌[16]和结肠癌[17]研究，验证了此算法的有效性。然而，目前鲜有研究将ESTIMATE算法应用于宫颈癌的免疫、基质评分及预后分析。鉴此，本研究通过癌症基因组图谱(The Cancer Genome Atlas,TCGA)数据库获得宫颈癌表达谱数据，通过ESTIMATE算法分析免疫、基质评分及其与宫颈癌患者生存预后的关系，并结合高通量基因表达汇编(Gene Expression Omnibus,GEO)数据库中的数据集(GSE63514),通过差异表达基因(differentially expressed genes,DEGs)分析、功能富集和生存分析，筛选出关键基因，并探讨关键基因在宫颈癌的表达情况及其与患者预后的关联性。从而为宫颈癌免疫相关的生物标志物的研究以及宫颈癌免疫治疗提供参考。现报道如下。

1 资料与方法

1.1数据获取和预处理 以“cervical cancer”“transcriptome profiling”“Gene Expression Quantification”为检索词，在TCGA数据库(https://tcga-data.nci.nih.gov/tcga/)下载信息数据集，包含306例宫颈癌组织和3例正常宫颈组织数据,同时下载与之对应的临床信息。将生存状态和生存时间缺失的患者信息删除，对收集到的数据集通过R语言软件“estimate”程序包(ESTIMATE算法，https://r-forge.r-project.org)进行免疫、基质评分，并应用R语言软件的“survival”程序包对免疫、基质进行生存分析。以“cervical cancer”作为关键词在GEO数据库(https://www.ncbi.nlm.nih.gov/geo)下载基因芯片数据集GSE63514，平台：Affymetrix Human Genome U133 Plus 2.0 Array，包含24例正常宫颈组织和28例宫颈癌组织数据信息。生存分析以基因表达水平的中位数作为界值分为高表达组和低表达组。

1.2DEGs分析 根据免疫、基质评分中位数将病例数据分为高、低免疫评分组和高、低基质评分组，应用R语言“limma”程序包对免疫和基质评分相关的基因做DEGs分析，获得上调及下调基因。通过R软件分析GSE63514数据集的芯片信息，将宫颈癌组织和正常宫颈组织数据进行比较以获得DEGs，并应用R软件制作火山图将其表示出来。通过韦恩图获取上调基因和下调基因的交集。DEGs的筛选条件：基因表达量变化倍数≥2,调整后P≤0.01。

1.3功能富集分析 应用R软件的“cluster profile”程序包对DEGs进行GO富集分析和KEGG通路分析，之后再进行基因集富集分析(gene set enrichment analysis,GSEA)。

1.4蛋白质互作(protein-protein interaction,PPI)网络分析 通过在线数据库STRING(http://string-db.org)对DEGs构建PPI网络，将最低关联度设置为0.850。通过R软件对分析结果进行可视化，最后选择前13位节点的DEGs作为关键基因，通过“survival”程序包对其进行生存分析。

1.5PPI网络中基因与免疫细胞的相关性分析 通过网站(http://Timer.cistrome.org/)下载免疫细胞和基因矩阵LM22文件，基于免疫细胞基因矩阵和TCGA数据基因矩阵，通过R软件运行CIBESORT源代码，获得样本免疫细胞矩阵，分析PPI网络中基因与免疫细胞的相关性，通过R软件程序包“ggplot2”绘制相关性热图。

1.6关键基因的表达情况及外部验证 通过UALCAN在线网站(http://ualcan.path.uab.edu/index)分析SELL、CXCL9基因在宫颈癌组织中的表达水平，并分析其与宫颈癌患者生存预后的关联性。通过GEPIA2(http://gepia.cancer-pku.c)和UALCAN(http://ualcan.path.uab.edu/index)在线工具对基因分子水平进行差异表达分析。根据TCGA数据库的基因矩阵，利用GSEA软件进行KEGG通路分析，筛选出CXCL9、SELL的富集信号通路。

2 结果

2.1宫颈癌免疫和基质评分与患者预后的关联性分析结果 基于TCGA数据库数据分析，高免疫评分组患者的生存情况优于低免疫评分组(log-rank检验：χ2=4.400，P=0.035，见图1ⓐ)。基质评分与患者的生存预后未显示出显著关联(log-rank检验：χ2=0.735，P=0.391，见图1ⓑ)。

ⓐ免疫评分；ⓓ基质评分

2.2DEGs分析结果 基于TCGA数据库数据，根据免疫评分高低进行筛选，共获得了1 063个DEGs，其中包括640个上调基因和423个下调基因(见图2ⓐ)。基于GEO数据库数据进行DEGs分析，共筛选出1 903个DEGs，包括了1 698个上调基因和205个下调基因(见图2ⓑ)。将两者上调、下调基因一起合并制作韦恩图，最终得到7个下调基因和100个上调基因(见图2ⓒⓓ)。

2.3富集分析结果 对韦恩图筛选出的107个DEGs进行GO富集分析，结果显示，DEGs在生物过程中参与免疫效应T细胞活性调节过程；在细胞组成中参与了质膜的外侧面和颗粒膜分泌等组成；在分子功能中参与细胞因子受体结合、受体配体活动和激活信号受体的活性等(见图3ⓐ)。KEGG通路分析和GSEA结果均提示，DEGs主要富集于细胞因子受体的相互作用信号通路和细胞黏附分子通路(见图3ⓑ)。

ⓐGO分析结果图；ⓑKEGG富集分析结果图

2.4PPI网络分析和生存分析结果 应用STRING在线工具构建了一个包括49个节点和134条边的PPI网络(见图4ⓐ)，并筛选出前30个关键基因(见图4ⓑ)。对前13个基因(节点数≥8)进行Kaplan-Meier生存分析，结果显示，只有SELL和CXCL9基因与宫颈癌患者的生存预后相关(log-rank检验：χ2=10.460，P=0.001；χ2=5.664，P=0.017)。见图4ⓒⓓ。

ⓐPPI网络图；ⓑ关联度排前30个基因的柱状图；ⓒ基于TCGA数据的CXCL9高、低表达组生存分析结果图；ⓓ基于TCGA数据的SELL高、低表达组生存分析结果图

2.5关键基因与免疫细胞的相关性分析结果 对以上筛选出的13个基因进行免疫细胞相关性分析，结果显示其与B细胞、T细胞和NK杀伤免疫细胞存在显著相关性(P<0.05)。见图5。

图5 关键基因与免疫细胞相关性分析结果图

2.6CXCL9和SELL基因的表达情况及外部验证结果 通过GEPIA2和UALCAN在线网站对CXCL9和SELL基因的表达情况进行分析，结果显示，宫颈癌组织中CXCL9呈高表达(P<0.05)，但SELL表达水平与正常组织或癌旁组织比较差异不显著(P>0.05)。见图6ⓐ～ⓓ。根据GEPIA2和UALCAN在线工具分析结果均显示SELL的表达情况与患者生存预后具有关联性，高表达组的预后更优(P<0.05)。但仅UALCAN的分析结果显示，CXCL9高表达组的生存预后优于低表达组(P<0.05)，GEPIA2未获得有统计学意义的结果(P>0.05)。见6ⓔ～ⓗ。

ⓐⓑ基于UALCAN在线网站分析CXCL9和SELL在宫颈癌组织中的表达水平；ⓒⓓ基于GEPIA2在线网站分析CXCL9和SELL在宫颈癌组织中的表达水平；ⓔⓕ基于UALCAN在线网站验证CXCL9、SELL表达水平与宫颈癌患者生存预后的关联性；ⓖⓗ基于GEPIA2在线网站验证CXCL9和SELL与宫颈癌患者生存预后的关联性

2.7CXCL9和SELL的GSEA结果 GSEA结果显示，CXCL9参与了免疫T细胞、免疫B细胞反应、趋化因子调控和NK杀伤细胞介导的细胞毒性信号通路；SELL参与了免疫T细胞、免疫B细胞反应和趋化因子调控等信号通路。

3 讨论

3.1本研究基于TCGA数据库宫颈癌患者的基因表达谱和临床信息，通过ESTIMATE算法获得免疫、基质评分，并分析了评分高低与宫颈癌患者的预后关系。根据免疫评分高低获得DEGs，与GEO数据库获得的DEGs取交集后得到7个下调基因和100个上调基因。对这些基因进行GO功能富集分析，发现其富集在免疫效应调节过程、T细胞活性、对病毒反应、白细胞黏附及T细胞活性调节等，均为免疫调节相关的生物过程。KEGG分析提示DEGs主要富集于细胞因子受体的相互作用信号通路和细胞黏附分子通路。本研究还应用生物信息学工具CIBESORT[18]分析DEGs与免疫细胞的相关性，发现这些基因与B细胞、T细胞和NK杀伤免疫细胞存在显著相关性。以上生物信息学分析结果提示宫颈癌DEGs参与免疫调节过程并与免疫细胞有显著相关性，免疫调节及免疫相关基因在宫颈癌的发生、发展过程中发挥了重要作用。

3.2本研究通过对这些DEGs进行PPI网络分析，筛选出30个关键基因，并发现CXCL9和SELL表达与宫颈癌患者生存相关。外部验证结果也显示，CXCL9基因在宫颈癌组织中高表达，且高表达者预后更优。SELL在宫颈癌组织中的表达未见明显升高，但高表达者具有较好的生存预后。GSEA结果显示，CXCL9参与了免疫T细胞、免疫B细胞反应、趋化因子调控和NK杀伤细胞介导的细胞毒性信号通路；SELL参与了免疫T细胞、免疫B细胞反应和趋化因子调控等信号通路。

3.3趋化因子是一种可溶性蛋白质，可分为4种亚型：CXC、CC、CX3C和C[19]，通过调节浓度梯度协调细胞迁移。在肿瘤发展的早期阶段，趋化因子塑造了TME的免疫景观，在宫颈癌的发生、发展中发挥了重要作用。有研究显示，CXCL3在宫颈癌中高表达，同时其可促进宫颈癌的恶性表型[20]。徐蕾蕾等[21]的研究发现，CXCL13高表达可以导致宫颈癌患者产生不良预后。CXCL12可以降低宫颈癌细胞对放射线反应敏感度，有望成为新的放射治疗靶点[22]。CXCL9是趋化因子CXC亚族的一员，也被称为γ-干扰素诱导的单因子，可由TME中的髓系细胞在炎症条件下产生，它吸引表达CXCR3受体的细胞，包括活化的T细胞和NK细胞，并已被证明在免疫检查点治疗的反应中发挥作用[23]。CXCL9表达可能与肿瘤组织中募集淋巴细胞至TME发挥抗肿瘤免疫应答功能密切相关[24]。CXCL9具有抗血管生成的作用，其可通过与内皮细胞作用，招募T细胞、NK细胞共同破坏肿瘤血管的形成[25]。有研究显示，CXCL9高表达于结直肠癌组织，且与临床病理相关，而CXCL9高表达者的5年无瘤生存期及总生存期明显长于低表达组[26]。高表达的CXCL9还与肾癌、胰腺导管腺癌等许多恶性肿瘤的总体生存期延长相关[27-28]。也有研究者通过对卵巢癌生物信息挖掘发现CXCL9 mRNA在卵巢癌组织中高表达，与卵巢癌患者较好的预后显著相关。有研究发现，在非小细胞肺癌中，CXCL9过表达可以降低肿瘤血管的形成进而抑制癌细胞的转移和生长[29]。CXCL9可以招募淋巴细胞发挥抗肿瘤免疫应答作用，亦可以抵抗血管形成，进而抑制肿瘤的远处转移和生长[24]。然而，截至目前，笔者在其他趋化因子与肿瘤的关系中尚未发现以上抗肿瘤作用，可能与我们尚未能对数据完全挖掘有关。

综上所述，本研究基于ESTIMATE算法，应用宫颈癌表达数据进行生物学信息学分析，发现免疫相关基因在宫颈癌的发生、发展过程中发挥重要作用，并筛选出与预后相关的免疫基因趋化因子CXCL9，为进一步研究提供基础，为宫颈癌的免疫治疗提供了潜在靶点。