基于生物信息學分析的非小細胞肺癌診斷預后相關基因的篩選

楊燕霞,金 蓮,王 欣,張 潔,柳小平*

(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的新策略。