基于生物信息學的結腸癌核心基因篩選及驗證

曾成， 祁國萍， 沈穎，2， 陸文斌，2， 鄧建忠，2， 劉遷，2， 金建華，2

(1. 江蘇大學附屬武進醫院腫瘤科，江蘇 常州 213017；2. 徐州醫科大學武進臨床學院腫瘤科，江蘇 常州 213017)

結直腸癌是常見的消化道惡性腫瘤，其全球發病率僅次于乳腺癌、肺癌，居第3位，死亡率居癌癥相關死亡的第2位[1]。目前認為結腸癌與遺傳、環境、飲食、炎癥性腸病、性別、種族等相關[2]，但其發病確切分子機制尚不清楚，臨床亦缺乏早期診斷和預后判斷的分子標志物。基于高通量測序及基因芯片技術的快速發展，有研究者通過生物信息學方法篩選結腸癌的潛在核心基因[3-4]，但由于數據集和篩選方法不同，篩選結果存在差異。本研究通過分析GEO數據庫中結腸癌數據集，篩選結腸癌核心基因，利用TCGA數據庫驗證核心基因表達水平并做生存分析，最后通過體外實驗探討預后相關基因對結腸癌細胞增殖的影響，旨在獲得更多與結腸癌發生發展相關的潛在核心基因。

1 材料和方法

1.1 數據集來源

從GEO數據庫(www.ncbi.nlm.nih.gov/geo/)中下載數據集，篩選條件：人類結腸癌全基因組表達譜數據，樣本中包含結腸癌組織及癌旁組織，樣本量>50。本研究選擇GSE23878、GSE37182和GSE74602數據集進行后續分析。GSE23878數據集基于GPL570平臺，包含結腸癌組織35例，癌旁組織24例；GSE37182數據集基于GPL6947平臺，包含結腸癌組織84例，癌旁組織88例；GSE74602數據集基于GPL6104平臺，包含結腸癌組織30例，癌旁組織30例。

1.2 細胞系及試劑

人結腸癌HCT116、SW620細胞購自中國科學院上海生物化學與細胞研究所。DMEM、RPMI 1640培養基、胎牛血清均為美國Gibco公司產品；pECMV、pECMV-AURKA、pEGFP、pEGFP-TIMP1質粒由淼靈質粒平臺構建；轉染試劑TurboFect Transfection、兔抗人Aurora A蛋白激酶(AURKA)一抗、兔抗人β-肌動蛋白一抗、山羊抗兔二抗、MTT試劑盒、5×上樣緩沖液、細胞裂解液均購自上海生工生物工程有限公司；兔抗人基質金屬蛋白酶抑制物1(TIMP1)一抗購自武漢三鷹生物技術有限公司；ECL化學發光液購自合肥蘭杰柯科技有限公司；BCA蛋白定量試劑盒購自上海碧云天生物技術有限公司。

1.3 篩選差異表達基因

在R軟件(版本：4.0.2)中利用ggplot2、limma、pheatmap等軟件包處理3個數據集。差異表達基因(differentially expressed genes，DEGs)篩選標準如下：錯誤發現率(false discovery rate，FDR)校正后P值<0.05，且|log2FC|>1。倍數變化(fold change，FC)表示DEGs的差值倍數。對3個數據集篩選出的DEGs取交集，利用VennDiagram軟件包作韋恩圖，得到共有DEGs。

1.4 共有DEGs的基因本體論富集分析和京都基因與基因組百科全書信號通路分析

在R軟件(版本：4.0.2)中利用clusterProfiler、org.Hs.eg.db、enrichplot、ggplot2等軟件包對共有DEGs進行基因本體論(gene ontology，GO)富集分析和京都基因與基因組百科全書(kyoto encyclopedia of genes and genomes，KEGG)信號通路分析，GO富集分析包括分子功能、生物學過程及細胞組成。

1.5 構建蛋白質相互作用網絡及篩選核心基因

利用STRING(https：//string-db.org/)在線數據庫對共有DEGs構建蛋白質相互作用(protein-protein interaction， PPI)網絡，設置有效結合分數>0.4。隨后將PPI網絡導入Cytoscape軟件(版本：3.7.2)，利用CytoHubba插件篩選核心基因。CytoHubba插件中有EPC、BottleNeck、EcCentricity、Closeness、Radiality、Betweenness、Stress、ClusteringCoefficient、MCC、DMNC、MNC和Degree 共12種分析方法[5]。每種分析方法選擇前10個基因，所有基因按照在12種分析方法中出現次數進行排序，最后選擇出現次數最多的10個基因作為核心基因，并構建核心基因的PPI網絡。

1.6 基于TCGA數據庫的核心基因驗證和生存分析

利用基于TCGA數據庫的GEPIA在線網站(http：//gepia.cancer-pku.cn/ index.html)驗證核心基因在結腸癌及癌旁組織表達水平，并對核心基因進行生存分析。

1.7 細胞培養與轉染

在37 ℃、5% CO2條件下，將HCT116、SW620細胞分別置于含10%胎牛血清DMEM、RPMI 1640培養基中培養。取對數生長期細胞接種于6孔板，待細胞融合度達80%時，按照轉染說明書，將pECMV(空載體)、pECMV-AURKA、pEGFP(空載體)、pEGFP-TIMP1質粒各2 μg分別轉染HCT116和SW620細胞。

1.8 蛋白質印跡法檢測細胞AURKA、TIMP1表達

收集“1.7”轉染48 h后各組細胞，加入含1% PMSF的RIPA裂解液，冰上裂解30 min；4 ℃、15 000 r/min離心20 min后取上清液，BCA法行蛋白定量；加入5×上樣緩沖液，100 ℃沸水浴5 min；分別取40 μg樣品行SDS-PAGE，80 V分離蛋白；200 mA 90 min冰上轉PVDF膜；5%脫脂牛奶室溫封閉2 h；PBST洗膜3次，5 min/次；加入相應的用5% 脫脂牛奶稀釋的一抗β-肌動蛋白(1 ∶4 000，內參)、AURKA和TIMP1(均為1 ∶1 000)， 4 ℃孵育過夜；PBST洗膜3次，5 min/次；加入用PBS稀釋的二抗(1 ∶5 000)室溫孵育1 h；PBST洗膜3次，5 min/次；用ECL化學發光液暗室曝光；用Image J 軟件處理分析蛋白條帶。

1.9 MTT法檢測細胞增殖

取“1.7”轉染后各組細胞懸液接種于96孔板，每孔2×103個，每組5個復孔。5% CO2、37 ℃分別培養24、48、72 h，每孔加入20 μL MTT(5 mg/mL)，37 ℃孵育4 h；棄上清液，每孔加入150 μL二甲基亞砜，置搖床上低速振蕩10 min，待結晶物充分溶解，于酶聯檢測儀490 nm波長處測定光密度(D)值，計算細胞增殖率。細胞增殖率(%)=(實驗組D值-空白組D值)/(對照組D值-空白組D值)×100%。

1.10 統計學方法

2 結果

2.1 DEGs篩選結果

通過R軟件分析，GSE23878數據集DEGs 1 411個，其中上調基因502個，下調基因909個；GSE37182數據集DEGs 627個，其中上調基因261個，下調基因366個；GSE74602數據集DEGs 1 485個，其中上調基因816個，下調基因669個。3個數據集DEGs通過VennDiagram軟件包取交集后得到270個共有DEGs(圖1)。

圖1 DEGs韋恩圖

2.2 共有DEGs的GO富集分析及KEGG信號通路分析

GO富集分析結果顯示，共有DEGs生物學過程主要與生長負調控、細胞外基質組織、細胞外結構組織等相關；細胞組成主要富集于含膠原的細胞外基質、收縮纖維、肌節等；分子功能主要富集于受體配體活性、信號轉導受體激活劑、糖胺聚糖結合等。KEGG信號通路富集結果表明，共有DEGs主要與礦物質吸收、Wnt、緊密連接、NF-κB、細胞周期、細胞黏附分子等信號通路相關。見圖2。

圖2 共有DEGs的GO富集分析及KEGG信號通路分析

2.3 PPI網絡構建及核心基因的篩選

利用STRING數據庫對共有DEGs進行PPI網絡分析，結果顯示，PPI網絡具有明顯的交互作用(P<1.0×10-16)。將結果導入Cytoscape軟件(圖3A)，利用插件CytoHubba中的12種算法篩選出重復次數最多的前10個核心基因，由于第11個基因重復出現的次數和第10個基因重復出現的次數相同，本研究也將其納入核心基因，故最終獲得11個核心基因。核心基因在12種算法中重復出現次數如下：MYC8次、基質金屬蛋白酶抑制物1(TIMP1)6次、泛素偶聯酶E2C(UBE2C)5次、小窩蛋白1(CAV1)5次、Y 染色體中性別決定區相關的高遷移率組框9(SOX9)5次、C-X-C型趨化因子配體12(CXCL12) 5次、Aurora A蛋白激酶(AURKA) 4次、Ⅰ型膠原蛋白α1鏈(COL1A1)4次、細胞周期分裂蛋白20(CDC20)4次、DNA拓撲異構酶Ⅱα(TOP2A)4次、著絲粒蛋白F(CENPF)4次。11個核心基因重新導入STRING數據庫后構建核心基因的PPI網絡(圖3B)，結果顯示有明顯交互作用(P=7.66×10-7)。

A:共有DEGs的PPI網絡；B:核心基因的PPI網絡

2.4 核心基因表達驗證及生存分析

GEPIA在線網站包含TCGA數據庫中275例結腸癌組織和41例癌旁組織相關信息。11個核心基因通過GEPIA網站驗證顯示，其中MYC、TIMP1、UBE2C、SOX9、AURKA、COL1A1、CDC20、TOP2A和CENPF共9個基因在結腸癌中呈高表達，而CAV1和CXCL12在結腸癌中呈低表達，差異有統計學意義(P<0.05)，見圖4。

圖4 核心基因在結腸癌和癌旁組織表達水平的驗證

此外，本研究也利用該網站對核心基因進行生存分析，發現TIMP1表達與結腸癌患者總體生存期(overall survival，OS)呈顯著負相關(P<0.05)，即TIMP1在結腸癌樣本中表達越高，結腸癌患者OS越短，AURKA表達與結腸癌患者OS呈顯著正相關(P<0.05)，而COL1A1表達與結腸癌患者無病生存期呈顯著負相關(P<0.05)，見圖5。

圖5 3個核心基因的生存分析

2.5 TIMP1或AURKA對結腸癌細胞增殖能力的影響

為進一步驗證影響結腸癌患者OS相關基因(TIMP1和AURKA)在結腸癌細胞中的功能，本研究分別將pECMV，pECMV-AURKA，pEGFP，pEGFP-TIMP1質粒轉染入HCT116和SW620細胞。蛋白質印跡結果顯示，HCT116、SW620細胞pECMV-AURKA組AURKA蛋白相對表達量顯著高于相應pECMV組(P均<0.01)；HCT116、SW620細胞 pEGFP-TIMP1組TIMP1蛋白表達相對量顯著高于相應pEGFP 組(P均<0.01)，見圖6。MTT結果顯示，HCT116、SW620細胞pECMV-AURKA 組72 h細胞增殖能力明顯高于pECMV組(t=4.039，5.731，P均<0.05)，24、48 h兩組間差異無統計學意義。HCT116細胞 pEGFP-TIMP1組48、72 h細胞增殖能力明顯高于pEGFP組(t=11.716，5.673，P均<0.01)，24 h兩組間差異無統計學意義。SW620細胞pEGFP-TIMP1組72 h細胞增殖能力明顯高于pEGFP組(t=5.920，P<0.01)，24、48 h兩組間差異無統計學意義(圖7)。

圖6 蛋白質印跡檢測HCT116和SW620細胞中AURKA和TIMP1蛋白表達

a：P<0.05，與同時間點pECMV組比較；b：P<0.01， 與同時間點pEGFP組比較

3 討論

本研究通過下載GEO數據庫中GSE23878、GSE37182和GSE74602數據集，經分析獲得270個共有DEGs。GO富集和KEGG信號通路分析顯示這些基因主要與生長負調控、受體配體活性、信號轉導受體激活劑、Wnt信號通路、NF-κB信號通路、細胞周期信號通路等有關。在篩選核心基因時，為了減少單一算法的局限性，本研究充分利用CytoHubba軟件中12種算法，將每一種算法的前10個核心基因按照重復出現的次數進行排序，選取重復次數最多的11個核心基因(MYC、TIMP1、UBE2C、CAV1、SOX9、CXCL12、AURKA、COL1A1、CDC20、TOP2A和CENPF)進行后續驗證分析。為了驗證GEO數據集篩選結果的準確性，將11個核心基因通過基于TCGA數據庫的GEPIA在線網站進行驗證，結果顯示11個核心基因在TCGA數據庫的表達水平與GEO數據庫篩選結果一致，其中MYC、TIMP1、UBE2C、SOX9、AURKA、COL1A1、CDC20、TOP2A和CENPF在結腸癌中高表達，而CAV1和CXCL12在結腸癌中低表達。

MYC蛋白是一種轉錄因子，在眾多腫瘤中發揮重要作用。UBE2C蛋白是泛素-蛋白酶體系統的重要組成部分，UBE2C高表達與多種腫瘤不良預后相關[6]，但其在結腸癌中的分子機制研究較少。CAV1蛋白是細胞膜主要支架蛋白，CAV1在結腸癌細胞中過表達能夠引起啟動子CpG位點低甲基化，進而促進細胞增殖[7]。SOX9蛋白是轉錄因子SOX家族成員之一，陳玉昌等[8]研究顯示，SOX9 mRNA和蛋白表達在結腸癌組織中均上調，且SOX9蛋白與腫瘤的分化程度、TNM分期和淋巴轉移有關。CXCL12蛋白是與CXCR4(G蛋白偶聯受體)特異性結合的一種趨化因子，CXCL12/CXCR4信號與多種腫瘤的侵襲、遷移密切相關[9-10]。TOP2A蛋白是一種酶蛋白，與細胞增殖、凋亡及有絲分裂相關。研究顯示，TOP2A敲減可通過影響凋亡蛋白(Bcl-2、Bax)、侵襲相關蛋白(MMP2、MMP9)表達抑制結腸癌細胞的增殖和侵襲[11]。CENPF蛋白是著絲粒蛋白家族成員之一，在有絲分裂及腫瘤中起重要作用，被認為是結腸腺癌的新型分子標志物[12]。CDC20蛋白是細胞周期相關蛋白，其異常表達使有絲分裂發生錯誤，從而導致一些癌基因的過表達及抑癌基因的失活，最終抑制腫瘤細胞凋亡和促進腫瘤細胞增殖、侵襲、轉移[13]。有研究顯示，CDC20高表達與結腸癌的臨床分期、病理分化程度和TNM分期有關，且CDC20高表達結腸癌患者OS較CDC20低表達患者短[14]。本研究顯示，CDC20高表達與結腸癌患者OS無統計學意義(P>0.05)，這可能與采用不同的生存分析方法有關。因此，CDC20對結腸癌患者預后的影響有待進一步研究。

COL1A1蛋白是膠原蛋白家族成員之一，是細胞外基質重要組成成分，COL1A1過表達可導致結腸癌細胞上皮-間質轉化，進而促進結腸癌細胞肝轉移[15]。本研究顯示COL1A1高表達與結腸癌患者無病生存期呈負相關，與丁志祥等[4]研究結果一致。AURKA蛋白在細胞有絲分裂過程中起重要作用，是調節細胞周期關鍵分子。研究顯示，使用AURKA蛋白抑制劑Alisertib可顯著抑制MYC驅動的結腸癌細胞增殖[16]，也能促進KRAS驅動的結腸癌細胞死亡[17]。本研究顯示AURKA過表達能夠促進結腸癌細胞增殖，但AURKA高表達結腸癌患者OS卻更長，這可能與AURKA高表達增加結腸癌細胞化療敏感性有關[18]。AURKA在結腸癌中呈高表達，可作為結腸癌的預后標志物[19]。TIMP1蛋白是組織金屬蛋白酶抑制劑-1，在抑制金屬蛋白酶介導的細胞外基質轉化方面發揮重要作用，進而參與腫瘤的侵襲和轉移。研究表明，TIMP1敲除后能夠通過FAK-PI3K/AKT 和 MAPK信號通路抑制結腸癌細胞的增殖和轉移，且與TNM分期、無病生存期、血管侵犯以及淋巴結轉移有關[20]。本研究通過生物信息學發現TIMP1表達量越高的結腸癌患者OS越短，MTT實驗顯示結腸癌細胞中過表達TIMP1能夠顯著促進細胞增殖，該結果與相關報道[20]一致，故認為TIMP1蛋白可能為潛在的結腸癌分子標志物。

綜上所述，本研究通過生物信息學方法篩選出可能的結腸癌核心基因，包括MYC、TIMP1、UBE2C、CAV1、SOX9、CXCL12、AURKA、COL1A1、CDC20、TOP2A和CENPF，其中AURKA和TIMP1表達變化與結腸癌患者預后相關，但仍需更多實驗研究進行驗證。