基于生物信息学分析的宫颈鳞状细胞癌预后相关基因的筛选

曲木金作

(四川省凉山州妇幼保健计划生育服务中心,中国四川凉山615000)

宫颈癌是全世界女性中第四大最常见的癌症,在发展中国家女性中死亡率很高[1]。作为最常见的宫颈癌组织类型,宫颈鳞状细胞癌(cervical squamous cell carcinoma,CESC)是严重威胁女性健康的疾病,每年造成约273 200人的死亡[2]。近年来,随着癌症筛查和各种治疗手段如手术、放疗和化疗的发展,CESC临床预后得到一定改善。然而,由于在疾病早期缺乏有效的诊断方法,CESC转移和复发的风险仍然很高,预后效果欠佳。越来越多的证据表明,多种基因的异常表达参与了CESC的发生和发展过程[3～5]。鉴于CESC的高发病率和高死亡率,早期发现和风险评估对改善CESC患者的预后显得尤为重要。因此,寻找新的诊断、预后和治疗靶点的生物标志物,以提高宫颈癌患者的生存率是必要和迫切的。

当前,生物信息学已被广泛用于筛选基因组水平的遗传改变[6～7]。在本研究中,我们基于基因表达谱芯片数据集,利用R软件筛选出CESC与正常组织之间的差异表达基因(differentially expressed genes,DEGs),并对其进行了功能和通路富集分析及蛋白质-蛋白质相互作用(protein-protein interaction,PPI)网络分析,随后对筛选出的hub基因进行了LASSO COX回归模型和总体生存率(overall survival,OS)分析,获得了7个与CESC患者生存密切相关的关键基因(ZWINT、CDC6、PBK、TOP2A、NUSAP1、CCNB1 和 CDK1)。我们的结果为CESC提供了潜在的预后生物标志物及治疗靶点,并为进一步研究CESC的分子机制提供了理论依据。

1 材料和方法

1.1 数据下载和预处理

通过GEO(Gene Expression Omnibus,http://www．ncbi．nlm．nih．gov/geo)和 TCGA(The Cancer Genome Atlas,https://cancergenome．nih．gov/)数据库下载CESC患者癌组织和癌旁组织的基因表达谱数据,其中GSE9750包括24例正常宫颈上皮组织样本和33例CESC样本,GSE63514包括24例正常宫颈组织样本和20例CESC样本,TCGA来源数据包含3例正常宫颈组织样本和306例CESC样本。所获取的数据集原始数据(CEL file)通过R语言(v3．6．1;http://r-project．org/)Affy 包[8]读取, 将原始的CEL文件去除批间差后进行背景校正、bootstrap校正、质量控制和标准化处理,并转化为探针表达矩阵。

1.2 差异表达基因的筛选

使用R语言limma包[9]筛选CESC样本和正常样本之间的DEGs。DEGs满足调整后的P＜0．05和|log2fold change(log2FC)|＞1。

1.3 DEGs的功能和通路富集分析

GO(Gene Ontology)是一种主要的生物信息学工具,用于注释基因并分析这些基因的生物学过程,其涵盖了生物学的3个方面:细胞组分(cellular component,CC)、分子功能(molecular function,MF)和生物过程(biological process,BP)[10]。为了分析DEGs的功能,我们使用R软件clusterProfiler包[11]对 DEGs进行 GO 功能富集分析,P＜0．05 被认为是显著富集,然后使用enrichplot包对富集结果进行可视化。

GSEA(Gene Set Enrichment Analysis)是使用预定义的基因集,将基因按照其在两类样本中的差异表达程度进行排序,然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集[12]。GSEA检测的是基因集合而不是单个基因的表达变化,因此可以包含细微的表达变化,进而得到更为理想的结果。我们采用GSEA(v6．3,http://software．broadinstitute．org/gsea/index．jsp)对通路进行富集分析,选择 c2．cp．kegg．v6．0．symbols．gmt 作为参考基因集,选择错误发现率（false discovery rate,FDR）＜0．25 且 P＜0．05 作为截止标准,结果使用 R语言进行美化及可视化。

1.4 PPI网络构建及hub基因的选择

STRING(http://www．string-db．org/)是评估蛋白质-蛋白质相互作用信息的系统生物学工具[13]。为了评估差异表达基因所编码蛋白质之间的互作关系,我们使用STRING数据库进行PPI分析。随后使用 Cytoscape 软件(v3．7．1)[14]将 PPI网络可视化。为了更好地提取有价值的信息,我们利用cyto-Hubba插件[15]来识别hub基因。通过cytoHubba插件选择最大相关标准中的前20个基因作为hub基因。

1.5 预后相关基因的筛选

为了筛选出与CESC预后强相关的基因,我们同时对hub基因进行单变量Cox回归分析及LASSO COX回归分析,其中单变量Cox回归根据基因表达量的中位值将基因分为高表达和低表达,筛选标准为P＜0．05,然后利用R包survival[16]以及survminer进行作图分析。

2 结果

2.1 CESC中DEGs的鉴定

对GSE9750和GSE63514两个数据集进行标准化处理,处理前后的对比结果以小提琴图的形式呈现(图1)。经过标准化处理后,两个数据集中各样本处于同一水平,表明其一致性较高。数据预处理后,我们利用R软件进行差异分析,GSE9750和GSE63514两个数据集分别获得522和985个DEGs,结果以火山图形式呈现(图2A,B);TCGA数据库来源的表达谱获得6 466个DEGs。将3个数据集的DEGs取交集,结果如图2C所示,最终获得167个DEGs。

图1 标准化处理前后样本表达量的小提琴图(A)GSE9750数据集标准化处理之前的小提琴图;(B)GSE9750数据集标准化处理之后的小提琴图;(C)GSE63514数据集标准化处理之前的小提琴图;(D)GSE63514数据集标准化处理之后的小提琴图。Fig.1 Violin diagram of the expression amount of sample before and after standardized treatment(A)Violin diagram before standardization of GSE9750 dataset;(B)Violin diagram after standardization of GSE9750 dataset;(C)Violin diagram before standardization of GSE63514 dataset;(D)Violin diagram after standardization of GSE63514 dataset．

图2 CESC组织与正常组织中的差异表达基因(A)GSE9750数据集中DEGs的火山图;(B)GSE63514数据集中DEGs的火山图;(C)GSE9750、GSE63514和TCGA 3个数据集中DEGs交集的韦恩图。Fig.2 DEGs between CESC and normal tissues(A)DEGs volcano map in GSE9750 dataset;(B)DEGs volcano map in GSE63514 dataset;(C)Intersection Venn diagrams of DEGs in GSE9750,GSE63514 and TCGA．

2.2 DEGs功能和通路富集分析

为了更深入地了解DEGs的生物学功能,我们运用R软件clusterProfiler包对DEGs进行GO功能富集分析。GO分析结果表明,DEGs的生物过程(BP)显著富集在染色体分离、核分裂、表皮细胞分化、DNA复制和姐妹染色单体分离;细胞组分(CC)主要涉及染色体区域、纺锤体、染色体着丝粒区域、桥粒和MCM(minichromosome maintenance)复合体;分子功能(MF)主要富集于染色质结合、G蛋白偶联受体结合、细胞因子受体结合、DNA依赖性ATP酶活性、蛋白酶结合和趋化因子活性(图3A)。进一步的GSEA分析结果显示,富集的通路主要涉及DNA复制和细胞周期(图3B),其中MCM家族在两个通路的信号转导过程中具有重要作用。除MCM家族之外,细胞周期通路中比较关键的分子还有CCNB1、CDC6和BUB1B。

2.3 PPI网络构建及hub基因的选择

使用STRING在线工具构建DEGs的PPI网络,并应用Cytoscape软件将其可视化,结果如图4A所示。此外,利用cytoHubba插件选择最大相关标准中的前20个基因作为hub基因,即MAD2L1、ZWINT、BUB1B、RRM2、TTK、AURKA、CDC6、PBK、RAD51AP1、DTL、TOP2A、KIF11、DLGAP5、KIF20A、NCAPG、RFC4、NUSAP1、CCNB1、MELK 及 CDK1为所得hub基因(图4B)。

2.4 预后相关基因的筛选

LASSO COX回归结果显示13个hub基因(MAD2L1、ZWINT、RRM2、TTK、CDC6、PBK、TOP-2A、KIF11、KIF20A、NCAPG、NUSAP1、CCNB1 和CDK1)与预后相关(图 5)。Kaplan-Meier曲线显示,ZWINT、DTL、CCNB1、CDC6、TOP2A、CDK1、PBK、RFC4及NUSAP1的低mRNA表达水平与CESC患者较差的生存预后相关(图6)。在上述结果中,ZWINT、CDC6、PBK、TOP2A、NUSAP1、CCNB1 及CDK1为LASSO COX回归和OS生存分析中重叠的hub基因,表明这7个hub基因可能是CESC异常信号传导途径中的关键参与者,可作为CESC潜在的预后生物标志物。

图3 差异表达基因的GO分析(A)及GSEA通路富集分析(B)Fig.3 GO analysis of DEGs(A)and enrichment analysis of GSEA pathway(B)

3 讨论

CESC是全世界女性中最常见的恶性肿瘤之一,具有较高的发病率[17]。CESC的复发和转移风险较高,早期诊断和治疗对于改善CESC患者的预后至关重要。因此,迫切需要探索可用于早期诊断、靶向治疗或预后评估的新型潜在生物标志物,以改善CESC患者预后。微阵列分析是一种高通量技术,可以同时检测数千个基因的表达水平。如今,基因的异常表达被认为是CESC发生和发展的因素之一,并且越来越多的研究表明CESC中一些失调的基因可能成为诊断和预后的候选生物标志物[4～5,18]。因此,我们通过生物信息学分析CESC的基因表达谱,探索其分子机制,鉴定可能作为CESC生物标志物和治疗靶点的重要分子。

图4 差异表达基因PPI网络和hub基因(A)PPI网络分析图。节点大小表示聚类系数,节点越大,聚类系数越大,基因在网络中占据的比重就越大;节点颜色表示度,度越大,该节点连线就越多,蓝色代表度较大,黄色居中,橙色最小;连线粗细表示综合得分,得分越高线越粗,越说明两个蛋白质之间存在互作;(B)Hub基因示意图。颜色越红,富集分数越高,颜色越黄,富集分数越小。Fig.4 PPI network of DEGs and hub genes(A)PPI network analysis graph．The node size represents the clustering coefficient．The larger the node size,the larger the clustering coefficient,and the greater the proportion of genes in the network．The node color represents the degree．The greater the degree,the more connections the node will have．Blue color means a greater degree,yellow means medium and orange smallest．The thickness of the connection lines represents the comprehensive score．The thicker the line,the higher the score,and the more likely the interaction between the two proteins will exist;(B)Hub gene schematic map．Darker red means the higher enrichment fraction,and yellow means relatively lower．

图5 20个hub基因的LASSO系数谱横轴越往左,自由度越小,γ越大,系数会越趋于0;不同颜色对应不同基因。Fig.5 LASSO coefficient profiles of the 20 hub genesThe degree of freedom towards further left of the transverse axis becomes smaller,and the larger the gamma,the closer the coefficient is to zero．Different colors correspond to different genes．

图6 TCGA队列中hub基因的Kaplan-Meier生存曲线黄线代表基因高表达组,绿线代表基因低表达组,P＜0．05具有统计学意义。Fig.6 Kaplan-Meier survival curves for hub genes in the TCGA cohortsThe yellow line represents high gene expression group,and the green line represents low gene expression group．P＜0．05 is statistically significant．

在本研究中,我们从GEO和TCGA数据库下载了CESC的基因表达谱数据集,并使用生物信息学方法对其进行了深入分析,获得了CESC组织和正常组织之间的DEGs。研究结果显示,在CESC组织和正常组织之间共鉴定出167个DEGs。GO功能富集分析显示,DEGs主要涉及染色体分离、核分裂、表皮细胞分化及DNA复制等生物过程,介导染色质结合、G蛋白偶联受体结合、细胞因子受体结合及DNA依赖性ATP酶活性等分子功能,同时这些基因的表达产物主要富集于染色体区域、纺锤体、染色体着丝粒区域、桥粒和MCM复合体。此外,GSEA通路富集结果显示,富集的通路主要涉及DNA复制和细胞周期。Zhu等[19]研究表明SNAP23通过诱导细胞周期G2/M期阻滞抑制宫颈癌的进展。另有研究报道,hnRNPA2/B1可通过抑制PI3K/Akt信号通路诱导细胞周期阻滞,进而抑制宫颈癌细胞增殖和侵袭[20]。MCM基因家族在DNA复制和细胞周期通路中具有重要作用,其编码的MCM蛋白家族是细胞复制周期中的重要调节因子,在判断肿瘤预后方面具有重要价值。例如:MCM2和MCM6的高表达与肝癌患者的预后不良相关[21];MCM5的高表达可能是非小细胞肺癌患者的独立不良预后生物标志物[22]。据报道,CCNB1、CDC6和BUB1B也与肿瘤的发生发展密切相关[23～25]。以上信息表明,我们的数据挖掘结果与已有研究结果相符。

此外,我们还构建了DEGs的PPI网络并筛选出了20个hub基因,分别为MAD2L1、ZWINT、BUB1B、RRM2、TTK、AURKA、CDC6、PBK、RAD51-AP1、DTL、TOP2A、KIF11、DLGAP5、KIF20A、NCAPG、RFC4、NUSAP1、CCNB1、MELK 和 CDK1。以上基因中的一部分已在先前的研究中被证明与CESC密切相关。例如:RRM2、AURKA和KIF20A的高表达与CESC较差的总生存率密切相关[26～28];TTK、BUB1B和MELK与宫颈癌的转移相关[25,29];KIF11在宫颈癌进展过程中介导胞质分裂[30];SIX1表达增加可导致RFC4的表达量上调,进而促进宫颈癌的发生、发展和侵袭性生长[31]。另外,Kim等[32]研究发现,MAD2L1在CESC中表达上调,表明其参与了宫颈癌的发生发展。然而,现阶段仍缺乏 RAD51AP1、DTL、DLGAP5 和 NCAPG 在CESC中的相关研究。

生物信息学背景下普遍存在着高维数据,所谓的“高维”即待估计的未知参数的个数是样本量的一个或几个数量级。以往大部分研究仅仅使用OS生存分析对疾病预后进行预测,而LASSO COX回归适用于分析高维度、强相关、小样本的生存资料数据,因此我们同时使用LASSO COX回归模型及OS生存分析来评估这20个hub基因对CESC患者存活的影响,以进一步提高预测结果的可信度。本研究结果显示,ZWINT、CDC6、PBK、TOP2A、NUSAP1、CCNB1 及 CDK1 的低mRNA表达水平与CESC患者较差的生存预后相关,表明这些基因可能在CESC的发生发展、侵袭转移及复发中起关键作用。Peres等[33]研究表明,TOP2A的表达增加可以促进宫颈癌细胞分裂,并且它可以用作宫颈癌免疫组织化学的生物标志物。CCNB1已被证明通过调节宫颈癌细胞的凋亡在宫颈癌的发生发展中发挥重要作用[23]。有研究报道,CDC6可能是宫颈高级别和侵袭性病变的生物标志物[24];PBK和NUSAP1的高表达与宫颈癌的转移相关[34～35]。Li等[36]报道 CDK1 是宫颈癌进展的关键效应因子,可能成为宫颈癌的潜在靶点。有意思的是,目前关于ZWINT对CESC的影响暂未见报道。对ZWINT进行深入的功能机制研究,可为阐明CESC的分子机制提供有价值的见解。

总之,我们使用生物信息学方法对CESC的3个基因表达谱数据集进行了深入挖掘,筛选出了20个hub基因,并通过进一步的LASSO COX回归及OS分析得到7个与CESC患者生存密切相关的关键基因(ZWINT、CDC6、PBK、TOP2A、NUSAP1、CCNB1和 CDK1),其中 CDC6、PBK、TOP2A、NUSAP1、CCNB1和CDK1已被证明与宫颈鳞状细胞癌密切相关,ZWINT在CESC的诊断、治疗靶点及预后方面的潜力有待进一步探索和验证。