杨燕霞,金 莲,王 欣,张 洁,柳小平*

(1.甘肃省第二人民医院呼吸科,中国甘肃兰州730030；2.重庆市北碚区中医院,中国重庆400700)

作为肺癌的主要类型,非小细胞肺癌(nonsmall cell lung cancer,NSCLC)约占肺癌的 85%[1]。尽管近年来在NSCLC的筛查、诊断及治疗等方面已经取得了一些进展,但患者的5年生存率仍然很低,只有17%[2]。随着基因组时代的开启及分子生物学研究的发展,在分子水平研究生命现象和疾病的发生机制引起了人们的极大关注,近年来NSCLC靶向新药开发的研究也越来越受到重视[3]。然而,多数患者会对靶向药物产生抗药性,导致靶向药物治疗效果欠佳。因此,发现新的NSCLC治疗靶点对提高NSCLC患者的生存率有着重要的意义。综合生物信息分析是将不同数据集整合到一起,获得更多的临床样本,以实现可信度更高的分析,从而为实验研究提供可行的思路[4～6]。本研究旨在利用生物信息学方法,从GEO及TCGA数据库中分别筛选与NSCLC相关的差异表达基因(differentially expressed genes,DEGs),然后取交集,对所获得的DEGs展开功能相关性分析,紧接着通过蛋白质网络互作选取hub基因,并对hub基因同时进行ROC诊断分析及LASSO生存分析,选取与预后密切相关的诊断标志物,以期为NSCLC提供预后生物标志物及治疗靶点,并为进一步研究NSCLC的分子机制提供新的思路。

1 材料与方法

1.1 数据下载与预处理

从 GEO 数据库(http：//www.ncbi.nlm.nih.gov/geo)[7]中筛选样本来源可靠的NSCLC表达谱数据集,使用 R 语言(version 3.6.1；http：//r-project.org/)GEOquery 包[8]下载并分析数据集 GSE18842[9]、GSE101929[10],两者均是基于GPL570平台的人肺组织表达谱数据集。其中,GSE18842数据集包括45例正常肺组织和46例NSCLC组织,GSE101929数据集包括34例正常肺组织和32例NSCLC组织。同样,使用R语言RTCGAToolbox包[11]在TCGA数据库(https：//www.cancer.gov/tcga)[12]下载 NSCLC mRNA基因表达数据以及临床数据,使用R语言软件整理相关表达及表型数据。

1.2 差异表达基因的筛选

通过affy包[13]将原始的CEL文件进行背景校正及均一化处理,校正效果使用密度图进行可视化,并转化为探针表达矩阵；根据Bioconductor平台对应的GPL平台注释文件,并对探针进行基因注释。通过limma软件包[14]筛选出DEGs,以|log2fold change(log2FC)|＞2且 P＜0.05为差异截取标准。为了展示DGEs的差异表达情况,以ggplot2软件包[15]绘制火山图。对于TCGA数据库中获取的数据,则使用edgeR包[16]进行DEGs的筛选,条件同样满足|log2FC|＞2且 P＜0.05。最后,使用 VennDiagram包[17]对三者取交集并以韦恩图进行可视化。

1.3 GSEA分析、差异表达基因的基因本体论和通路富集分析

基因集富集分析(gene set enrichment analysis,GSEA)[18]通过评估一个预先定义的基因集的基因在与表型相关度排序的基因表中的分布趋势,从而判断其对表型的贡献。基因本体论(gene ontology,GO)[19]是用来注释基因及其产物的常用方法,大规模基因的注释经常使用该分析方法。文中使用R语言clusterProfiler包[20]对DEGs进行GO和KEGG(kyoto encyclopedia of genes and ge-nomes)通路富集分析,同时以c2.cp.kegg.v6.0.symbols.gmt作为参考基因集进行GSEA分析,P＜0.05认为具有统计学意义。

1.4 蛋白质互作网络分析

STRING(version 11.0；http：//string-db.org)[21]是用于评估蛋白质-蛋白质相互作用(protein-protein interaction,PPI)信息的在线工具。Cytoscape常用于复杂网络的可视化[22],其插件cytoHubba可用于计算基因所得度值,常用来筛选hub基因[23]。首先将240个DEGs导入STRING中,得到它们的互作关系,再将所得互作关系导入Cytoscape软件,并用cytoHubba以Degree算法为标准筛选hub基因,定义得分排名前20的基因为所得hub基因。

1.5 Hub基因ROC诊断分析及LASSO生存分析

受试者操作特征曲线(receiver operating characteristic,ROC)能够直观地鉴别各诊断指标的诊断效能,ROC曲线越靠近左上角,曲线下面积(area under curve,AUC)越大,诊断价值越高[24]。文中使用R语言pROC包[25]对所获得的hub基因矩阵进行ROC诊断分析,以筛选具有诊断价值的hub基因；使用软件包glmnet[26]对hub基因进行LASSO COX回归分析,从hub基因中筛选出跟预后强相关的基因。

2 结果

2.1 差异表达基因筛选结果

GEO数据库来源数据集GSE18842、GSE1-01929的标准化处理结果如图1A所示,两组样本密度图曲线基本重合,可见两组样本来源可靠。GSE18842、GSE101929数据集DEGs火山图展示结果如图1B所示,GSE18842数据集中共筛选出735个DEGs,GSE101929数据集中共筛选出858个DEGs。此外,TCGA数据库中共筛选出951个DEGs。对GSE18842数据集、GSE101929数据集和TCGA数据库所得DEGs取交集,共筛选出240 个 DEGs(图 2)。

2.2 GSEA分析及GO和KEGG通路富集分析

GSEA分析不需要对基因进行表达差异的筛选,能保留表达变化不大但功能重要的基因,因此相比于GO和KEGG富集分析,该方法保留了更多的信息。本文的GSEA分析结果显示,NSCLC组富集的通路主要涉及DNA修复和细胞周期,其中MCM基因家族以及BUB基因家族在其中作用突出,结果如图3所示。进一步采用R语言对240个DEGs进行GO和KEGG通路富集分析,结果如图4所示。DEGs主要涉及核分裂、染色体分离的调控、核分裂调控等生物过程,主要与细胞外基质(extracellular matrix,ECM)受体相互作用信号通路以及细胞黏附分子、细胞周期等通路相关。

图1 数据集预处理及DEGs火山图(A)GSE18842数据集和GSE101929数据集标准化处理之后的密度图；(B)GSE18842数据集和GSE101929数据集中所获DEGs的火山图。Fig.1 Dataset preprocessing and volcano map of DEGs(A)The density map after standardization of GSE18842 and GSE101929 datasets；(B)The DEGs volcano map in GSE18842 and GSE101929 datasets.

图2 数据集中DEGs交集的Venn图Fig.2 Venn map of differential gene intersection of datasets

2.3 DEGs所编码蛋白质之间的相互作用分析

将240个DEGs输入STRING工具,然后将所得互作数据(图5A)导入Cytoscape中,使用插件cytoHubba找出hub基因,结果如图5B所示。IL6、UBE2C、KIAA0101、TOP2A、MAD2L1、CDC20、CEP55、RRM2、MKI67、CDC6、CCNB1、KIF2C、TPX2、AURKA、CDKN3、MELK、CDCA5、CENPF、NUF2、ZWINT为所得hub基因。

2.4 ROC诊断分析及LASSO生存分析

AUC＞0.5的情况下,AUC值越接近1,表明诊断标志物的诊断效果越好。基于GEO数据集,我们利用R语言绘制了20个hub基因的ROC曲线。结果如图6所示,hub基因的AUC基本位于 0.7～0.9,其中 UBE2C(AUC=0.939)、TOP2A(AUC=0.927)、RRM2(AUC=0.927)、CCNB1(AUC=0.928)、MKI67(AUC=0.930)、AURKA(AUC=0.931)和MELK(AUC=0.950)相对具有较高的诊断价值。利用LASSO回归分析hub基因对NSCLC预后的影响,结果如图 7所示,基因 IL6、KIAA0101、MKI67、TPX2、AURKA、CDKN3、CDCA5 均与 NSCLC患者的生存预后显著相关。

3 讨论

图3 GSEA通路富集分析Fig.3 Enrichment analysis of GSEA pathway

图4 差异表达基因的GO分析及KEGG通路富集分析Fig.4 GO analysis and KEGG pathway enrichment analysis of DEGs

NSCLC是导致全球癌症相关死亡的主要原因之一,其死亡率目前呈上升趋势[27]。虽然得益于免疫靶向治疗,不少患者的生活质量得到改善,但是晚期NSCLC患者的预后依旧很差。而且,尽管人们在NSCLC样本的基因组学、转录组学、蛋白质组学以及代谢组学等领域的研究中发现了不少有望成为诊断和预后的候选生物标志物[28～29],但是仍然没有改变当前NSCLC患者预后较差的事实。为提高肺癌患者的存活率,给更多的NSCLC患者带来希望,现阶段仍迫切需要发现新的有效的诊断和预后标记物。因此,进一步阐明NSCLC的发生发展机制,寻找有效的预后生物标志物,对NSCLC患者的预后意义重大。

本研究采用生物信息学的方法分析NSCLC和正常对照组之间的差异表达基因。分别从GEO、TCGA数据库下载NSCLC数据集,经预处理后,将数据集合并取交集,共获取240个DEGs。为了解NSCLC所涉及的通路,对DEGs进行了GSEA分析、GO功能注释和KEGG通路富集分析。GSEA结果显示DEGs主要富集到细胞周期及DNA修复通路,之前已有大量研究表明两者在NSCLC的发生发展及预后中起着重要的作用[30～34],这也进一步证实我们的数据整合分析结果是可靠的。GO分析结果显示,DEGs主要参与的生物过程为核分裂、染色体分离的调控、核分裂调控等；主要富集到细胞外基质；主要富集于细胞外基质结构成分以及胞外基质结合等。相关研究表明,核分裂、染色体分离的调控、核分裂调控与肿瘤的发生发展及转移有着密切关系[35～36]；细胞外基质结构成分以及胞外基质与肿瘤转移和侵袭有关[37]。此外,KEGG通路富集分析显示,ECM受体相互作用信号通路在NSCLC中起着一定作用,细胞与ECM之间的特异性相互作用由主要成分为整合素的跨膜分子介导,这些相互作用可以控制细胞黏附、迁移及周期[38～39]；而大量研究表明,整合素在NSCLC发生发展中起着重要的桥梁作用[40～42]。因此,ECM-受体相互作用信号通路有望成为NSCLC潜在的药物治疗靶点。由此可见,我们的研究结果与既往研究发现相一致。

图5 差异表达基因所编码蛋白质的PPI分析图(A)PPI网络分析图。节点的大小表示聚类系数,节点越大,聚类系数越大,说明该基因在网络中占据比重就越大。节点颜色表示度,度越大,说明该节点连线就越多,度从大到小分别为橙、黄、蓝。线的粗细代表综合得分,得分越高线越粗。线的颜色代表共表达,同一种颜色说明两蛋白间存在互作关系；(B)Hub基因所编码蛋白质的互作示意图。颜色越红越深,富集分数越高。Fig.5 Protein-protein interaction analysis of differential gene-encoded proteins(A)The PPI network analysis diagram.The size of a node represents the clustering coefficient,the larger the node,the larger the clustering coefficient,indicating that the gene occupies a larger proportion in the network.The node color indicates the degree,the greater the degree,the more connected the node.Degrees from big to small are orange,yellow,and blue.The thickness of the line represents the comprehensive score,the higher the score,the thicker the line.The color of the line represents co-expression,the same color indicates that there is an interaction between the two proteins；(B)The interaction diagram of hub geneencoded proteins.The darker the color is,the higher the concentration is.

图6 Hub基因ROC诊断曲线图横坐标为假阳性率,越接近零准确率越高；纵坐标为敏感度,也称为真阳性率,越大代表准确率越好。Fig.6 ROC diagnostic curves of hub genesThe abscissa shows the false positive rate；the closer it to zero,the higher the accuracy.The ordinate represents sensitivity,also known as the true positive rate.The higher the rate,the greater the accuracy.

图7 Hub基因LASSO回归图每一条曲线代表了一个自变量系数的变化轨迹,纵坐标是系数的值,下横坐标是log(lambda),上横坐标是此时模型中非零系数的个数。Fig.7 Hub gene LASSO regression mapEach curve represents the change track of the coefficient of each independent variable.The ordinate is the value of the coefficient.The lower abscissa is log(lambda),and the upper abscissa is the number of non-zero coefficients in the model at the time.

为进一步筛选与NSCLC诊断及预后密切相关的基因,我们利用STRING得到DEGs的互作网络关系,再利用Cytoscape插件cytoHubba分析得到了 20 个 hub 基因 IL6、UBE2C、KIAA0101、TOP2A、MAD2L1、CDC20、CEP55、RRM2、MKI67、CDC6、CCNB1、KIF2C、TPX2、AURKA、CDKN3、MELK、CDCA5、CENPF、NUF2、ZWINT。其中,IL6、UBE2C、KIAA0101、TOP2A、MAD2L1、RRM2、MKI67、CCNB1、TPX2、AURKA、CDKN3 在 NSCLC 发生发展中的作用研究较为深入,机制相对较为明确[43～53]。同时,CDC6、CEP55、MELK、CDCA5、CENPF、NUF2、CDC20在NSCLC的研究中也有见报道[54～57],但是,对于其具体导致NSCLC发生的分子机制仍有待进一步研究。ZWINT所编码蛋白质在间期细胞的细胞质中分布均匀,是动粒形成和纺锤体检查点活动所必需的[58]。2019年Peng等[59]的研究表明其可能是肺癌治疗的下一个重要靶点,因此,关于其参与NSCLC发生发展的机制研究值得进一步深入。KIF2C编码一种类似运动蛋白的蛋白质,该蛋白质作为一种依赖于微管的分子马达,能使正端微管解聚,从而促进有丝分裂染色体分离[60],而染色体分离与肿瘤的发生发展存在密切的关系,虽然目前尚无研究表明KIF2C与NSCLC有关系,但以上结果均提示KIF2C可能参与了NSCLC的发生发展。

此外,我们对20个hub基因进行了ROC诊断分析及LASSO生存分析,ROC诊断分析结果表明 UBE2C、TOP2A、RRM2、CCNB1、MKI67、AURKA、MELK相对具有较高的诊断价值,LASSO生存分析结果则显示 IL6、KIAA0101、MKI67、TPX2、AURKA、CDKN3、CDCA5 均与 NSCLC 患者生存预后显著相关,它们可能是NSCLC潜在的预后生物标志物。

综上所述,细胞周期和DNA修复对NSCLC发生发展起着关键作用；基于诊断分析、预后分析以及文献复习的结果,我们预测ZWINT、KIF2C、MELK、CDCA5可能在NSCLC中发挥着重要的作用；同时,ECM-受体相互作用信号通路与NSCLC密切相关,相关机制值得进一步深入研究。总之,这些结果为阐明NSCLC发生发展的分子机制提供了理论依据,并确定了ZWINT、KIF2C、MELK、CDCA5可能成为诊断生物标志物、潜在治疗靶点和预后指标的新关键基因,有助于开发诊断和治疗NSCLC的新策略。