基于生物信息學(xué)分析的宮頸鱗狀細(xì)胞癌預(yù)后相關(guān)基因的篩選

曲木金作

(四川省涼山州婦幼保健計(jì)劃生育服務(wù)中心,中國(guó)四川涼山615000)

宮頸癌是全世界女性中第四大最常見(jiàn)的癌癥,在發(fā)展中國(guó)家女性中死亡率很高[1]。作為最常見(jiàn)的宮頸癌組織類型,宮頸鱗狀細(xì)胞癌(cervical squamous cell carcinoma,CESC)是嚴(yán)重威脅女性健康的疾病,每年造成約273 200人的死亡[2]。近年來(lái),隨著癌癥篩查和各種治療手段如手術(shù)、放療和化療的發(fā)展,CESC臨床預(yù)后得到一定改善。然而,由于在疾病早期缺乏有效的診斷方法,CESC轉(zhuǎn)移和復(fù)發(fā)的風(fēng)險(xiǎn)仍然很高,預(yù)后效果欠佳。越來(lái)越多的證據(jù)表明,多種基因的異常表達(dá)參與了CESC的發(fā)生和發(fā)展過(guò)程[3～5]。鑒于CESC的高發(fā)病率和高死亡率,早期發(fā)現(xiàn)和風(fēng)險(xiǎn)評(píng)估對(duì)改善CESC患者的預(yù)后顯得尤為重要。因此,尋找新的診斷、預(yù)后和治療靶點(diǎn)的生物標(biāo)志物,以提高宮頸癌患者的生存率是必要和迫切的。

當(dāng)前,生物信息學(xué)已被廣泛用于篩選基因組水平的遺傳改變[6～7]。在本研究中,我們基于基因表達(dá)譜芯片數(shù)據(jù)集,利用R軟件篩選出CESC與正常組織之間的差異表達(dá)基因(differentially expressed genes,DEGs),并對(duì)其進(jìn)行了功能和通路富集分析及蛋白質(zhì)-蛋白質(zhì)相互作用(protein-protein interaction,PPI)網(wǎng)絡(luò)分析,隨后對(duì)篩選出的hub基因進(jìn)行了LASSO COX回歸模型和總體生存率(overall survival,OS)分析,獲得了7個(gè)與CESC患者生存密切相關(guān)的關(guān)鍵基因(ZWINT、CDC6、PBK、TOP2A、NUSAP1、CCNB1 和 CDK1)。我們的結(jié)果為CESC提供了潛在的預(yù)后生物標(biāo)志物及治療靶點(diǎn),并為進(jìn)一步研究CESC的分子機(jī)制提供了理論依據(jù)。

1 材料和方法

1.1 數(shù)據(jù)下載和預(yù)處理

通過(guò)GEO(Gene Expression Omnibus,http://www．ncbi．nlm．nih．gov/geo)和 TCGA(The Cancer Genome Atlas,https://cancergenome．nih．gov/)數(shù)據(jù)庫(kù)下載CESC患者癌組織和癌旁組織的基因表達(dá)譜數(shù)據(jù),其中GSE9750包括24例正常宮頸上皮組織樣本和33例CESC樣本,GSE63514包括24例正常宮頸組織樣本和20例CESC樣本,TCGA來(lái)源數(shù)據(jù)包含3例正常宮頸組織樣本和306例CESC樣本。所獲取的數(shù)據(jù)集原始數(shù)據(jù)(CEL file)通過(guò)R語(yǔ)言(v3．6．1;http://r-project．org/)Affy 包[8]讀取, 將原始的CEL文件去除批間差后進(jìn)行背景校正、bootstrap校正、質(zhì)量控制和標(biāo)準(zhǔn)化處理,并轉(zhuǎn)化為探針表達(dá)矩陣。

1.2 差異表達(dá)基因的篩選

使用R語(yǔ)言limma包[9]篩選CESC樣本和正常樣本之間的DEGs。DEGs滿足調(diào)整后的P＜0．05和|log2fold change(log2FC)|＞1。

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

GO(Gene Ontology)是一種主要的生物信息學(xué)工具,用于注釋基因并分析這些基因的生物學(xué)過(guò)程,其涵蓋了生物學(xué)的3個(gè)方面:細(xì)胞組分(cellular component,CC)、分子功能(molecular function,MF)和生物過(guò)程(biological process,BP)[10]。為了分析DEGs的功能,我們使用R軟件clusterProfiler包[11]對(duì) DEGs進(jìn)行 GO 功能富集分析,P＜0．05 被認(rèn)為是顯著富集,然后使用enrichplot包對(duì)富集結(jié)果進(jìn)行可視化。

GSEA(Gene Set Enrichment Analysis)是使用預(yù)定義的基因集,將基因按照其在兩類樣本中的差異表達(dá)程度進(jìn)行排序,然后檢驗(yàn)預(yù)先設(shè)定的基因集合是否在這個(gè)排序表的頂端或者底端富集[12]。GSEA檢測(cè)的是基因集合而不是單個(gè)基因的表達(dá)變化,因此可以包含細(xì)微的表達(dá)變化,進(jìn)而得到更為理想的結(jié)果。我們采用GSEA(v6．3,http://software．broadinstitute．org/gsea/index．jsp)對(duì)通路進(jìn)行富集分析,選擇 c2．cp．kegg．v6．0．symbols．gmt 作為參考基因集,選擇錯(cuò)誤發(fā)現(xiàn)率（false discovery rate,FDR）＜0．25 且 P＜0．05 作為截止標(biāo)準(zhǔn),結(jié)果使用 R語(yǔ)言進(jìn)行美化及可視化。

1.4 PPI網(wǎng)絡(luò)構(gòu)建及hub基因的選擇

STRING(http://www．string-db．org/)是評(píng)估蛋白質(zhì)-蛋白質(zhì)相互作用信息的系統(tǒng)生物學(xué)工具[13]。為了評(píng)估差異表達(dá)基因所編碼蛋白質(zhì)之間的互作關(guān)系,我們使用STRING數(shù)據(jù)庫(kù)進(jìn)行PPI分析。隨后使用 Cytoscape 軟件(v3．7．1)[14]將 PPI網(wǎng)絡(luò)可視化。為了更好地提取有價(jià)值的信息,我們利用cyto-Hubba插件[15]來(lái)識(shí)別hub基因。通過(guò)cytoHubba插件選擇最大相關(guān)標(biāo)準(zhǔn)中的前20個(gè)基因作為hub基因。

1.5 預(yù)后相關(guān)基因的篩選

為了篩選出與CESC預(yù)后強(qiáng)相關(guān)的基因,我們同時(shí)對(duì)hub基因進(jìn)行單變量Cox回歸分析及LASSO COX回歸分析,其中單變量Cox回歸根據(jù)基因表達(dá)量的中位值將基因分為高表達(dá)和低表達(dá),篩選標(biāo)準(zhǔn)為P＜0．05,然后利用R包survival[16]以及survminer進(jìn)行作圖分析。

2 結(jié)果

2.1 CESC中DEGs的鑒定

對(duì)GSE9750和GSE63514兩個(gè)數(shù)據(jù)集進(jìn)行標(biāo)準(zhǔn)化處理,處理前后的對(duì)比結(jié)果以小提琴圖的形式呈現(xiàn)(圖1)。經(jīng)過(guò)標(biāo)準(zhǔn)化處理后,兩個(gè)數(shù)據(jù)集中各樣本處于同一水平,表明其一致性較高。數(shù)據(jù)預(yù)處理后,我們利用R軟件進(jìn)行差異分析,GSE9750和GSE63514兩個(gè)數(shù)據(jù)集分別獲得522和985個(gè)DEGs,結(jié)果以火山圖形式呈現(xiàn)(圖2A,B);TCGA數(shù)據(jù)庫(kù)來(lái)源的表達(dá)譜獲得6 466個(gè)DEGs。將3個(gè)數(shù)據(jù)集的DEGs取交集,結(jié)果如圖2C所示,最終獲得167個(gè)DEGs。

圖1 標(biāo)準(zhǔn)化處理前后樣本表達(dá)量的小提琴圖(A)GSE9750數(shù)據(jù)集標(biāo)準(zhǔn)化處理之前的小提琴圖;(B)GSE9750數(shù)據(jù)集標(biāo)準(zhǔn)化處理之后的小提琴圖;(C)GSE63514數(shù)據(jù)集標(biāo)準(zhǔn)化處理之前的小提琴圖;(D)GSE63514數(shù)據(jù)集標(biāo)準(zhǔn)化處理之后的小提琴圖。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組織與正常組織中的差異表達(dá)基因(A)GSE9750數(shù)據(jù)集中DEGs的火山圖;(B)GSE63514數(shù)據(jù)集中DEGs的火山圖;(C)GSE9750、GSE63514和TCGA 3個(gè)數(shù)據(jù)集中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的生物學(xué)功能,我們運(yùn)用R軟件clusterProfiler包對(duì)DEGs進(jìn)行GO功能富集分析。GO分析結(jié)果表明,DEGs的生物過(guò)程(BP)顯著富集在染色體分離、核分裂、表皮細(xì)胞分化、DNA復(fù)制和姐妹染色單體分離;細(xì)胞組分(CC)主要涉及染色體區(qū)域、紡錘體、染色體著絲粒區(qū)域、橋粒和MCM(minichromosome maintenance)復(fù)合體;分子功能(MF)主要富集于染色質(zhì)結(jié)合、G蛋白偶聯(lián)受體結(jié)合、細(xì)胞因子受體結(jié)合、DNA依賴性ATP酶活性、蛋白酶結(jié)合和趨化因子活性(圖3A)。進(jìn)一步的GSEA分析結(jié)果顯示,富集的通路主要涉及DNA復(fù)制和細(xì)胞周期(圖3B),其中MCM家族在兩個(gè)通路的信號(hào)轉(zhuǎn)導(dǎo)過(guò)程中具有重要作用。除MCM家族之外,細(xì)胞周期通路中比較關(guān)鍵的分子還有CCNB1、CDC6和BUB1B。

2.3 PPI網(wǎng)絡(luò)構(gòu)建及hub基因的選擇

使用STRING在線工具構(gòu)建DEGs的PPI網(wǎng)絡(luò),并應(yīng)用Cytoscape軟件將其可視化,結(jié)果如圖4A所示。此外,利用cytoHubba插件選擇最大相關(guān)標(biāo)準(zhǔn)中的前20個(gè)基因作為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 預(yù)后相關(guān)基因的篩選

LASSO COX回歸結(jié)果顯示13個(gè)hub基因(MAD2L1、ZWINT、RRM2、TTK、CDC6、PBK、TOP-2A、KIF11、KIF20A、NCAPG、NUSAP1、CCNB1 和CDK1)與預(yù)后相關(guān)(圖 5)。Kaplan-Meier曲線顯示,ZWINT、DTL、CCNB1、CDC6、TOP2A、CDK1、PBK、RFC4及NUSAP1的低mRNA表達(dá)水平與CESC患者較差的生存預(yù)后相關(guān)(圖6)。在上述結(jié)果中,ZWINT、CDC6、PBK、TOP2A、NUSAP1、CCNB1 及CDK1為L(zhǎng)ASSO COX回歸和OS生存分析中重疊的hub基因,表明這7個(gè)hub基因可能是CESC異常信號(hào)傳導(dǎo)途徑中的關(guān)鍵參與者,可作為CESC潛在的預(yù)后生物標(biāo)志物。

圖3 差異表達(dá)基因的GO分析(A)及GSEA通路富集分析(B)Fig.3 GO analysis of DEGs(A)and enrichment analysis of GSEA pathway(B)

3 討論

CESC是全世界女性中最常見(jiàn)的惡性腫瘤之一,具有較高的發(fā)病率[17]。CESC的復(fù)發(fā)和轉(zhuǎn)移風(fēng)險(xiǎn)較高,早期診斷和治療對(duì)于改善CESC患者的預(yù)后至關(guān)重要。因此,迫切需要探索可用于早期診斷、靶向治療或預(yù)后評(píng)估的新型潛在生物標(biāo)志物,以改善CESC患者預(yù)后。微陣列分析是一種高通量技術(shù),可以同時(shí)檢測(cè)數(shù)千個(gè)基因的表達(dá)水平。如今,基因的異常表達(dá)被認(rèn)為是CESC發(fā)生和發(fā)展的因素之一,并且越來(lái)越多的研究表明CESC中一些失調(diào)的基因可能成為診斷和預(yù)后的候選生物標(biāo)志物[4～5,18]。因此,我們通過(guò)生物信息學(xué)分析CESC的基因表達(dá)譜,探索其分子機(jī)制,鑒定可能作為CESC生物標(biāo)志物和治療靶點(diǎn)的重要分子。

圖4 差異表達(dá)基因PPI網(wǎng)絡(luò)和hub基因(A)PPI網(wǎng)絡(luò)分析圖。節(jié)點(diǎn)大小表示聚類系數(shù),節(jié)點(diǎn)越大,聚類系數(shù)越大,基因在網(wǎng)絡(luò)中占據(jù)的比重就越大;節(jié)點(diǎn)顏色表示度,度越大,該節(jié)點(diǎn)連線就越多,藍(lán)色代表度較大,黃色居中,橙色最小;連線粗細(xì)表示綜合得分,得分越高線越粗,越說(shuō)明兩個(gè)蛋白質(zhì)之間存在互作;(B)Hub基因示意圖。顏色越紅,富集分?jǐn)?shù)越高,顏色越黃,富集分?jǐn)?shù)越小。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個(gè)hub基因的LASSO系數(shù)譜橫軸越往左,自由度越小,γ越大,系數(shù)會(huì)越趨于0;不同顏色對(duì)應(yīng)不同基因。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隊(duì)列中hub基因的Kaplan-Meier生存曲線黃線代表基因高表達(dá)組,綠線代表基因低表達(dá)組,P＜0．05具有統(tǒng)計(jì)學(xué)意義。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數(shù)據(jù)庫(kù)下載了CESC的基因表達(dá)譜數(shù)據(jù)集,并使用生物信息學(xué)方法對(duì)其進(jìn)行了深入分析,獲得了CESC組織和正常組織之間的DEGs。研究結(jié)果顯示,在CESC組織和正常組織之間共鑒定出167個(gè)DEGs。GO功能富集分析顯示,DEGs主要涉及染色體分離、核分裂、表皮細(xì)胞分化及DNA復(fù)制等生物過(guò)程,介導(dǎo)染色質(zhì)結(jié)合、G蛋白偶聯(lián)受體結(jié)合、細(xì)胞因子受體結(jié)合及DNA依賴性ATP酶活性等分子功能,同時(shí)這些基因的表達(dá)產(chǎn)物主要富集于染色體區(qū)域、紡錘體、染色體著絲粒區(qū)域、橋粒和MCM復(fù)合體。此外,GSEA通路富集結(jié)果顯示,富集的通路主要涉及DNA復(fù)制和細(xì)胞周期。Zhu等[19]研究表明SNAP23通過(guò)誘導(dǎo)細(xì)胞周期G2/M期阻滯抑制宮頸癌的進(jìn)展。另有研究報(bào)道,hnRNPA2/B1可通過(guò)抑制PI3K/Akt信號(hào)通路誘導(dǎo)細(xì)胞周期阻滯,進(jìn)而抑制宮頸癌細(xì)胞增殖和侵襲[20]。MCM基因家族在DNA復(fù)制和細(xì)胞周期通路中具有重要作用,其編碼的MCM蛋白家族是細(xì)胞復(fù)制周期中的重要調(diào)節(jié)因子,在判斷腫瘤預(yù)后方面具有重要價(jià)值。例如:MCM2和MCM6的高表達(dá)與肝癌患者的預(yù)后不良相關(guān)[21];MCM5的高表達(dá)可能是非小細(xì)胞肺癌患者的獨(dú)立不良預(yù)后生物標(biāo)志物[22]。據(jù)報(bào)道,CCNB1、CDC6和BUB1B也與腫瘤的發(fā)生發(fā)展密切相關(guān)[23～25]。以上信息表明,我們的數(shù)據(jù)挖掘結(jié)果與已有研究結(jié)果相符。

此外,我們還構(gòu)建了DEGs的PPI網(wǎng)絡(luò)并篩選出了20個(gè)hub基因,分別為MAD2L1、ZWINT、BUB1B、RRM2、TTK、AURKA、CDC6、PBK、RAD51-AP1、DTL、TOP2A、KIF11、DLGAP5、KIF20A、NCAPG、RFC4、NUSAP1、CCNB1、MELK 和 CDK1。以上基因中的一部分已在先前的研究中被證明與CESC密切相關(guān)。例如:RRM2、AURKA和KIF20A的高表達(dá)與CESC較差的總生存率密切相關(guān)[26～28];TTK、BUB1B和MELK與宮頸癌的轉(zhuǎn)移相關(guān)[25,29];KIF11在宮頸癌進(jìn)展過(guò)程中介導(dǎo)胞質(zhì)分裂[30];SIX1表達(dá)增加可導(dǎo)致RFC4的表達(dá)量上調(diào),進(jìn)而促進(jìn)宮頸癌的發(fā)生、發(fā)展和侵襲性生長(zhǎng)[31]。另外,Kim等[32]研究發(fā)現(xiàn),MAD2L1在CESC中表達(dá)上調(diào),表明其參與了宮頸癌的發(fā)生發(fā)展。然而,現(xiàn)階段仍缺乏 RAD51AP1、DTL、DLGAP5 和 NCAPG 在CESC中的相關(guān)研究。

生物信息學(xué)背景下普遍存在著高維數(shù)據(jù),所謂的“高維”即待估計(jì)的未知參數(shù)的個(gè)數(shù)是樣本量的一個(gè)或幾個(gè)數(shù)量級(jí)。以往大部分研究?jī)H僅使用OS生存分析對(duì)疾病預(yù)后進(jìn)行預(yù)測(cè),而LASSO COX回歸適用于分析高維度、強(qiáng)相關(guān)、小樣本的生存資料數(shù)據(jù),因此我們同時(shí)使用LASSO COX回歸模型及OS生存分析來(lái)評(píng)估這20個(gè)hub基因?qū)ESC患者存活的影響,以進(jìn)一步提高預(yù)測(cè)結(jié)果的可信度。本研究結(jié)果顯示,ZWINT、CDC6、PBK、TOP2A、NUSAP1、CCNB1 及 CDK1 的低mRNA表達(dá)水平與CESC患者較差的生存預(yù)后相關(guān),表明這些基因可能在CESC的發(fā)生發(fā)展、侵襲轉(zhuǎn)移及復(fù)發(fā)中起關(guān)鍵作用。Peres等[33]研究表明,TOP2A的表達(dá)增加可以促進(jìn)宮頸癌細(xì)胞分裂,并且它可以用作宮頸癌免疫組織化學(xué)的生物標(biāo)志物。CCNB1已被證明通過(guò)調(diào)節(jié)宮頸癌細(xì)胞的凋亡在宮頸癌的發(fā)生發(fā)展中發(fā)揮重要作用[23]。有研究報(bào)道,CDC6可能是宮頸高級(jí)別和侵襲性病變的生物標(biāo)志物[24];PBK和NUSAP1的高表達(dá)與宮頸癌的轉(zhuǎn)移相關(guān)[34～35]。Li等[36]報(bào)道 CDK1 是宮頸癌進(jìn)展的關(guān)鍵效應(yīng)因子,可能成為宮頸癌的潛在靶點(diǎn)。有意思的是,目前關(guān)于ZWINT對(duì)CESC的影響暫未見(jiàn)報(bào)道。對(duì)ZWINT進(jìn)行深入的功能機(jī)制研究,可為闡明CESC的分子機(jī)制提供有價(jià)值的見(jiàn)解。

總之,我們使用生物信息學(xué)方法對(duì)CESC的3個(gè)基因表達(dá)譜數(shù)據(jù)集進(jìn)行了深入挖掘,篩選出了20個(gè)hub基因,并通過(guò)進(jìn)一步的LASSO COX回歸及OS分析得到7個(gè)與CESC患者生存密切相關(guān)的關(guān)鍵基因(ZWINT、CDC6、PBK、TOP2A、NUSAP1、CCNB1和 CDK1),其中 CDC6、PBK、TOP2A、NUSAP1、CCNB1和CDK1已被證明與宮頸鱗狀細(xì)胞癌密切相關(guān),ZWINT在CESC的診斷、治療靶點(diǎn)及預(yù)后方面的潛力有待進(jìn)一步探索和驗(yàn)證。