生物信息學分析篩選結直腸癌靶基因及評估預后價值

車運誠，陳 梅，張 昱,2，張文靜,3*

(1．昆明理工大學 醫(yī)學院， 昆明 650500；2．云南省第一人民醫(yī)院 消化內(nèi)科， 昆明 650032；3．云南省第一人民醫(yī)院 腫瘤內(nèi)科， 昆明 650032)

結直腸癌(Colorectal cancer，CRC)是最常見的消化道惡性腫瘤，也是癌癥患者死亡的主要原因之一，嚴重危害著人類健康與生存。據(jù)報道，2018年全球結直腸癌新發(fā)病例180萬例，但死亡病例高達88.1萬例，僅次于肺癌[1-2]。而且，約15%～25%的結直腸癌患者在確診時即合并多發(fā)肝轉移，其中絕大多數(shù)患者的肝轉移灶無法獲得根治性切除。目前，手術和全身化療仍是結直腸癌的主要治療方法[3]。近年來，以抗表皮生長因子受體(Epidermal growth factor receptor EGFR，如西妥昔單抗、帕尼單抗)和抗血管內(nèi)皮生長因子(Vascular endothelial growth factor VEGF，如貝伐單抗，其他小分子抑制劑如瑞戈非尼等)為代表的靶向藥物以及程序性死亡受體1/程序性死亡-配體1(programmed cell death-1/ programmed cell death-Ligand 1, PD-1/PD-L1)免疫治療的應用，為晚期結直腸癌患者帶來了希望[4-5]。但上述精準治療都存在優(yōu)勢人群的篩選、顯著獲益的人數(shù)比例低，而且都存在一定的耐藥性。因此，尋找結直腸癌發(fā)生與發(fā)展過程的關鍵基因和通路，有助于我們認識結直腸癌潛在的發(fā)病機制，為其診治和靶點治療提供新的參考。

腫瘤精準治療的前提在于獲得腫瘤分子圖譜和分類，從中篩選并獲得臨床治療靶點；而后者得益于二代測序技術的進步和分子靶向藥物的問世。基于快速發(fā)展的微陣列和高通量測序技術，生物信息學的發(fā)展也日益蓬勃。公共數(shù)據(jù)庫如GEO和TCGA數(shù)據(jù)庫存儲大量、復雜的生物信息數(shù)據(jù)，通過對這些數(shù)據(jù)進行差異分析、聚類比對、生物分子網(wǎng)絡和通路分析以及可視化作圖等能夠更好地理解腫瘤的發(fā)病機制。

利用生物信息學的初步挖掘，篩選出潛在的結直腸癌發(fā)生發(fā)展相關功能基因，希望以此獲得更多與結直腸癌發(fā)生、進展相關的分子機制，從而為挖掘結直腸潛在的預后指標和腫瘤靶點提供一定的幫助。

1 材料和方法

1.1 材料

在NCBI(National center for biotechnology)的公共GEO數(shù)據(jù)庫(http://www.ncbi.nlm.nih.gov/geo/, Gene expression omnibus)中檢索含有人源結直腸癌樣本信息的數(shù)據(jù)芯片，下載數(shù)據(jù)集GSE106582[6]作為分析對象。芯片信息：Illumina HumanHT-12 V4.0 expression beadchip，平臺是GPL10558，該芯片中包含77例結直腸癌組織和117例癌旁組織(距離癌組織5cm處的粘膜組織)的表達數(shù)據(jù)。

1.2 篩選差異表達基因

通過R包對樣本進行PCA主成分分析，將樣本分為癌組織和癌旁組織，利用GEO數(shù)據(jù)庫的分析工具GEO2R對GSE106582芯片進行數(shù)據(jù)下載及差異表達基因篩選。篩選的標準為調(diào)準后p<0.05，且|log2FC|>1.5，并將探針名稱轉化為標準基因名稱。上調(diào)差異基因和下調(diào)差異基因結果以癌組織VS癌旁組織體現(xiàn)。

1.3 差異表達基因的GO功能富集和KEGG通路富集分析

使用在線分析工具DAVID 6.8(http://david.ncifcrf.gov/)對差異表達基因進行GO功能富集分析和KEGG通路富集分析。P<0.05認為具有統(tǒng)計學意義。

1.4 差異表達基因的相互作用分析

將差異表達基因導入在線數(shù)據(jù)庫STRING 10.5(https://string-db.org/)，構建結直腸癌差異基因蛋白質相互作用網(wǎng)絡，并運用Cytoscape3.1.1軟件進行可視化分析。使用軟件中的插入式分子復合物檢測(MCODE，使用默認參數(shù))篩選出 PPI 網(wǎng)絡中的樞紐模塊和關鍵基因。

1.5 差異表達基因的預后價值分析及ROC診斷

利用在線工具GEPIA2(http://gepia2.cancer-pku.cn/)分析關鍵基因表達水平與結直腸癌患者總生存期(Overall survival，OS)的相關性。根據(jù)結直腸癌患者差異表達基因的表達情況及中位值分為高表達組和低表達組，繪制生存曲線，若 log-rankP<0.05，則該基因被視為具有預后價值的基因。利用R語言包(PROC)對具有預后價值的基因進行受試者工作特征曲線(Receiver operating characteristic curve，ROC)分析，AUC>0.9則認為該基因對預后具有較高的靈敏度和特異度。下載數(shù)據(jù)集GSE21510對實驗結果進行驗證。

2 結果分析

2.1 樣本分組

對數(shù)據(jù)集GSE106582的樣本進行分組，其PCA結果顯示癌組織與癌旁組織的兩類樣本是相互獨立的，可將樣本分為癌組織和癌旁組織進行后續(xù)分析(見圖1)。

圖1 GSE106582 PCA分析圖Fig.1 PCA analysis diagram of GSE106582

2.2 數(shù)據(jù)處理與差異基因分析

在|log2FC|>1.5且調(diào)整后P<0.05的條件下，從數(shù)據(jù)集GSE106582中共得到199個差異表達基因，其中56個是上調(diào)基因，143個是下調(diào)基因(見圖2)。

圖2 基因的表達水平及分布火山圖Fig.2 Volcano map of the expression level and distribution of genes

2.3 GO富集分析

運用DAVID在線工具對差異表達基因進行GO富集分析發(fā)現(xiàn)，其上調(diào)基因主要與膠原蛋白分解代謝過程(Collagen catabolic process)，細胞外基質分解(Extracellular matrix disassembly)，細胞外基質組織(Extracellular matrix organization)和膠原纖維組織(Collagen fibril organization)等生物學過程有關(見圖3a)；而下調(diào)基因主要與碳酸氫鹽運輸(Bicarbonate transport)，負生長調(diào)控(Negative regulation of growth)，細胞對鋅離子的反應(Cellular response to zinc ion)和一碳代謝過程(One-carbon metabolic process)等生物學過程有關(見圖3b)。

圖3 差異表達基因GO富集過程結果Fig.3 Results of GO enrichment process of differentially expressed genes

2.4 KEGG富集分析

對差異表達基因進行KEGG通路富集分析發(fā)現(xiàn)，其上調(diào)基因主要參與細胞外基質-受體相互作用(ECM-receptor interaction)，蛋白質消化與吸收(Protein digestion and absorption)，PI3K/AKT信號通路(PI3K-Akt signaling pathway)和黏著斑(Focal adhesion)等信號通路(見圖4a)；而下調(diào)基因主要參與礦物質吸收(Mineral absorption)，藥物代謝-細胞色素P450(Drug metabolism-cytochrome P450)，化學致癌途徑(Chemical carcinogenesis)和氮代謝(Nitrogen metabolism)等信號通路(見圖4b)。

圖4 差異表達基因通路注釋結果Fig.4 Differentially expressed gene pathway annotation results

2.5 差異表達基因蛋白相互作用網(wǎng)絡的構建與分析

將差異表達基因導入在線分析數(shù)據(jù)庫STRING，并設定置信度 0.4 作為判斷相互作用是否有意義的標準，去除網(wǎng)絡中無連接的節(jié)點，構建PPI 網(wǎng)絡。將STRING構建的PPI網(wǎng)絡導入Cytoscape3.1.1軟件進行可視化分析，該網(wǎng)絡有152個節(jié)點，377條邊組成(圖5a)。使用MCODE插件分析出最顯著的相互作用模塊，該模塊網(wǎng)絡由13個節(jié)點，63條邊組成(見圖5b)。上述13個節(jié)點/基因分別為BGN，COL1A1，COL1A2，MMP3，MMP1，TIMP1，CXCL12，PLAU，CXCL1，CXCL8，CXCL5，MMP7，SPP1。其中1個為下調(diào)基因(CXCL12)，其余12個均為上調(diào)基因，后者有可能在結直腸的發(fā)生、發(fā)展過程中扮演癌基因角色，從而與結直腸癌的發(fā)生、發(fā)展與預后相關。

圖5 差異表達基因蛋白相互作用網(wǎng)絡構建及可視化結果Fig.5 Construction and visualization results of differentially expressed gene protein interaction network

2.6 預后價值分析

利用在線工具GEPIA2對顯著模塊中13個差異表達基因進行預后價值分析，結果顯示：其中4個基因BGN，COL1A2，TIMP1，SPP1的高表達組的結直腸癌患者的OS顯著低于相應的低表達組(見圖6a)。對上述的4個基因利用R語言包進行ROC診斷，結果顯示，4個基因BGN，COL1A2，TIMP1，SPP1的ROC曲線下面積AUC分別為0.928，0.925，0.931和0.888，除了SPP1基因之外，另外3個基因的AUC均大于0.9，表明這3個基因BGN，COL1A2，TIMP1可作為靈敏度和特異度較高的結直腸癌不良預后指標(見圖6b)。通過GSE21510進行驗證，BGN、COL1A2和TIMP1的診斷靈敏度要高于SPP1，與GSE106582實驗結果是一致的。

圖6 預后價值分析結果Fig.6 Prognostic value analysis results

3 討 論

通過生物信息學分析方法對GEO數(shù)據(jù)庫中的GSE106582 芯片篩選結直腸癌與癌旁組織的差異表達基因，共獲得199個差異基因，其中54個為上調(diào)基因，143個為下調(diào)基因基因。通過構建蛋白質相互作用網(wǎng)絡圖，得到13個核心差異表達基因，再應用STRING和GEPIA2工具對顯著模塊中13個差異表達基因進行預后價值分析，篩選出BGN、COL1A2和TIMP1可能是結直腸癌的不良預后指標。

雙糖鏈蛋白聚糖(biglycan，BGN)是富含亮氨酸的蛋白多糖家族中的關鍵成員之一，通常在神經(jīng)、骨、軟骨、皮膚和肌肉中表達，發(fā)揮調(diào)節(jié)上皮細胞的形態(tài)、生長、粘附、骨礦化、炎癥、遷移和分化功能[7]。BGN在膀胱癌、結直腸癌、胃癌等多種癌癥中表達水平較高[8-10]，參與多個腫瘤生物學過程包括細胞增殖、侵襲、上皮-間質轉化、血管生成和化療抵抗[9,11]。其表達水平受多個因素多個信號通路調(diào)控，包括p38[12]、HIF-1[13]和NF-κB[14]等。在結直腸癌中，BGN/TLR4/NF-κB軸在炎癥刺激下，能夠以正反饋方式最終增強NF-κB信號通路活性，從而促進癌癥的發(fā)生和進展[14]。

組織金屬蛋白酶抑制劑-1(Tissue inhibitor of metalloproteinases-1，TIMP1)是基質金屬蛋白酶家族(MMPs)的內(nèi)源性抑制劑。TIMP1的過表達與多種癌癥的不良預后相關，包括結直腸癌[15]、胃癌[16]、肺癌[17]。Batra等人發(fā)現(xiàn)，TIMP1的過度表達能夠增加增殖和凋亡相關基因的表達水平[18]；另有報道指出，通過誘導TIMP1特異性調(diào)控FAK-PI3K/AKT和MAPK通路，能夠抑制TIMP1的表達，從而降低腫瘤細胞的增殖和轉移，同時促進細胞凋亡[19]；高水平的癌胚抗原和TIMP1與較短的結直腸癌患者生存時間相關，尤其是血漿TIMP1表達水平可以作為結直腸癌的潛在不良預后指標[20]。

Ⅰ型膠原α2鏈(Collagen type1 alpha2，COL1A2) 屬于膠原蛋白基因家族中的成員之一。作為細胞外基質的主要成分之一，COL1A2主要通過細胞外基質受體通路和局部黏附通路影響細胞的增殖、分化、黏附和轉移[21]。現(xiàn)有研究發(fā)現(xiàn)，COL1A2在肝癌、卵巢癌等癌組織中的表達水平明顯高于正常組織，且高表達組患者的總體生存期較短[22-23]。在結直腸癌中，COL1A2的過度表達促進了EMT的發(fā)生。Mori[24]等報道 COL1A2 通過CpG島過甲基化引起膀胱癌細胞增殖和遷移。類似地，COL1A2過甲基化現(xiàn)象也存在于結直腸癌，且COL1A2的過表達能夠通過抑制NF-κB信號通路從而調(diào)控結直腸癌細胞的增殖、遷移與轉移[25]。以上報道均提示BGN、COL1A2和TIMP1是結直腸癌的不良預后指標，與本研究生物信息學預測分析結論一致。

通過GO和KEGG富集分析發(fā)現(xiàn)，上調(diào)的差異表達基因主要富集在與膠原蛋白分解代謝過程，細胞外基質分解，細胞外基質受體相互作用和PI3K/AKT信號通路。細胞外基質在調(diào)節(jié)細胞內(nèi)穩(wěn)態(tài)和細胞間相互作用中發(fā)揮一定的作用，同時膠原蛋白是ECM的主要組分，在惡性腫瘤的細胞粘附、形態(tài)改變、侵襲和遷移過程中扮演重要的角色[26]。PI3K/Akt信號通路是多種腫瘤發(fā)生發(fā)展和炎癥反應的重要信號通路，亦在促進腫瘤細胞增殖、調(diào)節(jié)細胞代謝、腫瘤生長、遷移和細胞骨架重塑等方面發(fā)揮重要作用[27]。下調(diào)的差異表達基因主要富集在碳酸氫鹽運輸、一碳代謝代謝過程、礦物質吸收、藥物代謝-細胞色素P450和氮代謝通路等生物過程。一方面，上述途徑主要集中在營養(yǎng)物質吸收和外源物質代謝等；另一方面，腸上皮細胞是防止細菌內(nèi)毒素等有害物質進入人體的重要屏障，其主要功能是吸收營養(yǎng)、毒素和藥物[28]。最近的研究表明，腸上皮細胞在維持腸道免疫穩(wěn)態(tài)方面起著重要的作用，而相關信號通路的異常表達已被報道與結直腸癌的發(fā)生有關[29]。分析結果為PI3K/Akt信號通路在結直腸癌中的重要作用提供了更多的理論依據(jù)。目前PI3K/Akt/mTOR通路抑制劑依維莫司已被批準臨床應用于晚期乳腺癌，而在結直腸癌中尚無報道。

4 結 論

綜上所述，本研究通過多種生物信息學分析方法確定了199個差異表達基因，并發(fā)現(xiàn)BGN、COL1A2和TIMP1在結直腸癌中表達上調(diào)，可能是結直腸癌獨立的不良預后因素。PI3K/Akt信號通路參與了結直腸癌的發(fā)生和轉移。然而，BGN、COL1A2和TIMP1的表達水平或其生物學功能機制是否與PI3K/Akt信號通路有關目前尚不明確。對進一步深入研究結直腸癌發(fā)生和進展的分子機制、預后指標及治療靶點的篩選具有重要意義，但需要進一步的實驗和臨床大樣本的驗證。