基于生物信息學分析的宮頸鱗狀細胞癌預后相關基因的篩選

曲木金作

(四川省涼山州婦幼保健計劃生育服務中心,中國四川涼山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的診斷、治療靶點及預后方面的潛力有待進一步探索和驗證。