李卓陽(yáng)， 張皓旻， 劉格良， 陳浩然， 智 鵬， 陳熙勐， 盧學(xué)春， 賀培鳳

1.山西醫(yī)科大學(xué)管理學(xué)院，山西 太原 030001；2.中國(guó)人民解放軍總醫(yī)院第二醫(yī)學(xué)中心血液科 國(guó)家老年疾病臨床醫(yī)學(xué)研究中心

結(jié)腸癌是最常見的惡性腫瘤之一，是我國(guó)惡性腫瘤死亡的第4常見原因[1]。結(jié)腸腺癌(colon adenocarcinoma，COAD)是發(fā)生于腺上皮細(xì)胞的惡性腫瘤，是結(jié)腸癌最主要的病理類型之一，目前的治療手段包括外科手術(shù)切除、化學(xué)療法、免疫療法等。由于復(fù)發(fā)等因素，術(shù)后患者的5年生存率僅為60%～70%[1-3]，提高患者的生存優(yōu)勢(shì)仍是一個(gè)挑戰(zhàn)。

全轉(zhuǎn)錄組測(cè)序(RNA-seq)可同時(shí)檢測(cè)全部已知基因的表達(dá)模式，通過(guò)差異表達(dá)基因分析，可鑒定腫瘤與正常組織之間的差異表達(dá)基因[4]。近年來(lái)，研究者基于癌癥基因圖譜(The Cancer Genome Atlas，TCGA)、基因表達(dá)綜合(Gene Expression Omnibus，GEO)等開源數(shù)據(jù)庫(kù)開展的生物信息分析發(fā)現(xiàn)，COAD的發(fā)生與大量基因的顯著變化有關(guān)。同時(shí)，基因間調(diào)控機(jī)制復(fù)雜，且其發(fā)展涉及多種信號(hào)通路異常[5-7]。目前，結(jié)腸癌已有預(yù)測(cè)預(yù)后基因的篩選研究，但未見基于基因表達(dá)與臨床表型關(guān)系篩選預(yù)后基因的相關(guān)研究。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis，WGCNA)是一種根據(jù)表達(dá)譜數(shù)據(jù)探索特定基因模塊與臨床表型間相關(guān)關(guān)系的系統(tǒng)生物學(xué)方法，可用于鑒定候選生物標(biāo)志物、預(yù)后基因及治療靶標(biāo)[8-9]。

本研究應(yīng)用WGCNA方法，對(duì)來(lái)自TCGA和GEO數(shù)據(jù)庫(kù)的COAD RNA-seq數(shù)據(jù)，構(gòu)建COAD的差異基因共表達(dá)網(wǎng)絡(luò)，篩選與COAD發(fā)生發(fā)展密切相關(guān)的基因模塊，隨后結(jié)合Kaplan-Meier方法鑒定COAD預(yù)后基因，為COAD預(yù)后相關(guān)基因的進(jìn)一步基礎(chǔ)及臨床研究提供方向和指導(dǎo)。

1 材料與方法

1.1 數(shù)據(jù)獲取從TCGA數(shù)據(jù)庫(kù)[10]中獲取COAD轉(zhuǎn)錄組數(shù)據(jù)和相應(yīng)的臨床信息。其中，轉(zhuǎn)錄組數(shù)據(jù)包括398例COAD樣本和39例正常結(jié)腸組織樣本；臨床信息包括COAD患者的生存狀態(tài)和生存時(shí)間。

從GEO數(shù)據(jù)庫(kù)[11]中以“colon adenocarcinoma”檢索COAD基因表達(dá)譜數(shù)據(jù)。篩選條件包括：研究物種為人或小鼠；全基因組表達(dá)芯片數(shù)據(jù)或轉(zhuǎn)錄組測(cè)序數(shù)據(jù)；生物學(xué)樣本組織來(lái)源類型一致；有對(duì)照組；每組至少3個(gè)生物學(xué)重復(fù)；實(shí)驗(yàn)設(shè)計(jì)思路清晰以及數(shù)據(jù)質(zhì)量良好。最終篩選得到GSE110224數(shù)據(jù)集[12]。該數(shù)據(jù)集包括17例原發(fā)性COAD樣本和17例配對(duì)的正常結(jié)腸組織樣本，基于GPL570平臺(tái)分析原發(fā)性COAD的整體基因表達(dá)變化。

1.2 差異表達(dá)基因分析采用R語(yǔ)言軟件包limma，分別對(duì)來(lái)自TCGA和GEO的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行數(shù)據(jù)標(biāo)準(zhǔn)化和差異表達(dá)基因(differentially expressed genes，DEGs)分析，得到兩組DEGs。差異表達(dá)基因的篩選標(biāo)準(zhǔn)，TCGA為|logFC|≥1、FDR<0.05；GEO為|logFC|≥1、FDR<0.05。當(dāng)多個(gè)探針與一個(gè)相同的基因匹配時(shí)，則以平均值作為該基因的表達(dá)值。最后采用R語(yǔ)言軟件包ggplot2對(duì)DEGs的表達(dá)模式進(jìn)行可視化。

1.3 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析采用R語(yǔ)言軟件包WGCNA，分別對(duì)來(lái)自TCGA和GEO的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析。填補(bǔ)缺失值后，構(gòu)建鄰接矩陣并轉(zhuǎn)換為拓?fù)渲丿B矩陣。采用動(dòng)態(tài)剪切法構(gòu)建基因聚類樹狀圖，將表達(dá)相似的基因聚類為不同的基因共表達(dá)模塊，規(guī)定每個(gè)模塊最少基因數(shù)目為50。隨后采用Pearson相關(guān)分析計(jì)算每個(gè)基因與各個(gè)模塊、不同臨床特征基因顯著性的相關(guān)系數(shù)r值和P值，獲取腫瘤相關(guān)性最高的模塊及該模塊中的基因。最后采用R語(yǔ)言軟件包VennDiagram，將相關(guān)系數(shù)最高模塊中的基因與TCGA、GEO的DEGs取交集，得到最終的DEGs，用以后續(xù)分析。

1.4 GO富集和KEGG通路富集分析采用R語(yǔ)言軟件包c(diǎn)lusterProfiler，對(duì)上述DEGs進(jìn)行GO[13-14]富集分析和KEGG通路[15]富集分析。以FDR<0.05、P<0.05為閾值，篩選富集的GO條目和KEGG通路。

1.5 蛋白互作網(wǎng)絡(luò)構(gòu)建及COAD預(yù)后相關(guān)的核心基因篩選采用STRING v11.0開源數(shù)據(jù)庫(kù)(https://string-db.org/cgi/input.pl)構(gòu)建蛋白互作(protein-protein interaction，PPI)網(wǎng)絡(luò)，以描述DEGs編碼蛋白之間的相互作用關(guān)系。選擇物種為“Homosapiens”，設(shè)置可靠性閾值>0.4，去除游離節(jié)點(diǎn)后下載PPI網(wǎng)絡(luò)數(shù)據(jù)。將數(shù)據(jù)導(dǎo)入Cytoscape軟件對(duì)PPI網(wǎng)絡(luò)進(jìn)行可視化，并使用Cytohubba插件中的MMC算法，篩選值最高的10個(gè)節(jié)點(diǎn)作為COAD預(yù)后相關(guān)的核心基因。

1.6 生存分析確定COAD預(yù)后相關(guān)的關(guān)鍵基因采用R軟件語(yǔ)言包survival，基于TCGA數(shù)據(jù)庫(kù)的COAD基因表達(dá)譜數(shù)據(jù)和臨床信息(患者生存狀態(tài)和生存時(shí)間)，通過(guò)Kaplan-Meier生存分析法，分析與患者總生存期(overall survival, OS)顯著相關(guān)的基因，以初步篩選出與患者預(yù)后不良有關(guān)的基因。隨后采用在線分析工具GEPIA2(http://gepia2.cancer-pku.cn/)，分析上述10個(gè)核心基因與患者無(wú)病生存期(disease-free survival，DFS)之間的關(guān)系。以P<0.05為閾值，篩選同時(shí)與OS和DFS有關(guān)的基因，作為COAD預(yù)后相關(guān)的關(guān)鍵基因。

1.7 HPA數(shù)據(jù)庫(kù)驗(yàn)證關(guān)鍵基因的蛋白質(zhì)表達(dá)人類蛋白質(zhì)圖譜(The Human Protein Atlas，HPA)(https://www.proteinatlas.org/)數(shù)據(jù)庫(kù)是利用轉(zhuǎn)錄組學(xué)和蛋白質(zhì)組學(xué)技術(shù)，從RNA和蛋白水平研究人類不同組織和器官中的蛋白表達(dá)情況。采用HPA數(shù)據(jù)庫(kù)，分析關(guān)鍵基因所編碼的蛋白在腫瘤組織和正常組織中的不同表達(dá)之處，進(jìn)一步明確與COAD預(yù)后相關(guān)的關(guān)鍵基因。

2 結(jié)果

2.1 差異表達(dá)基因篩選從TCGA數(shù)據(jù)庫(kù)下載COAD數(shù)據(jù)集，經(jīng)分析得到差異表達(dá)基因3 544個(gè)，其中上調(diào)基因1 293個(gè)，下調(diào)基因2 251個(gè)。從GEO數(shù)據(jù)庫(kù)下載GSE110224數(shù)據(jù)集，經(jīng)分析得到差異表達(dá)基因515個(gè)，其中上調(diào)基因223個(gè)，下調(diào)基因292個(gè)(見圖1)。

圖1 TCGA和GEO中COAD的DEGs火山圖

2.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析基于加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析，來(lái)自TCGA和GEO的所有基因分別被分為20個(gè)和18個(gè)模塊(見圖2～3)。如圖所示，與腫瘤正相關(guān)性最強(qiáng)的模塊分別是TCGA黃色(r=0.54，P=1e-34)和GEO棕色(r=0.65，P=3e-05)，與腫瘤負(fù)相關(guān)性最強(qiáng)的模塊分別是TCGA棕色(r=-0.87，P=1e-135)和GEO青色(r=-0.63，P=7e-05)。分別獲取模塊TCGA黃色和GEO棕色、TCGA棕色和GEO青色的共同基因，將兩組基因合并后再與TCGA和GEO的DEGs取交集，得到最終的153個(gè)DEGs，進(jìn)行后續(xù)分析(見圖4)。

注：A：基因聚類樹狀圖；B：基因模塊與腫瘤之間相關(guān)性的熱圖。

注：A：基因聚類樹狀圖；B：特征基因模塊與COAD關(guān)系圖。

圖4 DEGs與WGCNA模塊中基因的韋恩圖

2.3 GO富集和KEGG通路富集分析GO富集從三個(gè)方面注釋了基因的生物學(xué)特性：生物學(xué)進(jìn)程(biological process，BP)、細(xì)胞成分(molecular function，CC)和分子功能(molecular function，MF)(見圖5)。其中，BP集中于離子運(yùn)輸、激素代謝等；CC主要與細(xì)胞的頂端部分、頂質(zhì)膜、刷狀緣、微絨毛等細(xì)胞結(jié)構(gòu)有關(guān)；MF則主要與離子和分子跨膜轉(zhuǎn)運(yùn)蛋白，以及碳酸鹽脫水酶活性、類固醇脫水酶活性等有關(guān)。

注：橫坐標(biāo)代表富集在各GO條目的基因數(shù)目比率，縱坐標(biāo)代表GO條目名稱。

篩選得到KEGG通路富集有21條(見圖6)，發(fā)現(xiàn)這些基因主要參與膽汁分泌、類固醇激素合成、戊糖和葡萄糖醛酸酯相互轉(zhuǎn)化，以及氮、視黃醇、絡(luò)氨酸、丙酮酸等代謝過(guò)程。

注：橫坐標(biāo)代表富集在各通路的基因數(shù)目比率，縱坐標(biāo)代表通路名稱。

2.4 PPI網(wǎng)絡(luò)構(gòu)建和COAD預(yù)后相關(guān)核心基因篩選PPI網(wǎng)絡(luò)中共出現(xiàn)153個(gè)節(jié)點(diǎn)和385條連線(見圖7A)。根據(jù)Closeness算法，凝聚素Ⅰ復(fù)合物亞基G(non-SMC condensin I complex subunit G，NCAPG)、細(xì)胞分裂周期6(cell division cycle 6 homolog，CDC6)、核受體亞家族1，組H，成員4(nuclear receptor subfamily 1, group H, member 4，NR1H4)、甲狀腺素受體結(jié)合因子13(thyroid hormone receptor interactor 13，TRIP13)、氯離子通道輔助蛋白1(chloride channel, calcium activated, family member 1，CLCA1)、腸促胰高素樣肽1類似物胰高血糖素(glucagon，GCG)、核苷酸還原酶M2肽(ribonucleotide reductase M2 polypeptide，RRM2)、蔗糖酶異麥芽糖酶復(fù)合物(sucrase isomaltase，SI)、周期蛋白B1(Cyclin B1，CCNB1)、叉頭框M1(forkhead box M1，F(xiàn)OXM1)等10個(gè)基因?yàn)槠渲械暮诵幕?見圖7B)。

注：A：153個(gè)DEGs的PPI網(wǎng)絡(luò)；B：153個(gè)DEGs的核心基因。

2.5 關(guān)鍵基因的確定與驗(yàn)證10個(gè)核心基因中，OS分析結(jié)果顯示，CLCA1低表達(dá)患者的OS明顯比高表達(dá)組短(P<0.001，見圖8)；DFS分析結(jié)果顯示，TRIP3低表達(dá)、CLCA1低表達(dá)與患者的DFS顯著相關(guān)(P<0.05，見圖9)。選擇CLCA1基因作為關(guān)鍵基因。根據(jù)HPA數(shù)據(jù)庫(kù)，與正常結(jié)腸組織相比，腫瘤組織中CLCA1基因的蛋白質(zhì)水平顯著降低(見圖10)。

圖8 OS分析結(jié)果

圖9 10個(gè)核心基因的DFS分析結(jié)果

注：A;正常結(jié)腸組織，B;結(jié)腸腺癌組織。

3 討論

本研究通過(guò)對(duì)TCGA和GEO數(shù)據(jù)庫(kù)中的COAD轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行差異表達(dá)基因分析和加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析，發(fā)現(xiàn)COAD主要與患者體內(nèi)的153個(gè)基因異常表達(dá)有關(guān)。富集分析發(fā)現(xiàn)這些基因大多與離子運(yùn)輸、激素代謝等生命活動(dòng)有關(guān)。最后通過(guò)蛋白互作網(wǎng)絡(luò)分析和生存分析，發(fā)現(xiàn)CLCA1基因與COAD患者的不良預(yù)后顯著相關(guān)。

細(xì)胞內(nèi)離子通道在所有細(xì)胞中無(wú)處不在，研究表明包括K+、Cl-、Ca2+和Na+在內(nèi)的離子通道在胃腸道癌癥中均有表達(dá)和失調(diào)，這可能是導(dǎo)致正常細(xì)胞向癌細(xì)胞轉(zhuǎn)化的重要原因[16]。本研究中DEGs主要與無(wú)機(jī)和有機(jī)離子通道的激活和抑制等生物學(xué)進(jìn)程有關(guān)，提示離子通道的異常表達(dá)或功能障礙對(duì)COAD患者體內(nèi)癌細(xì)胞轉(zhuǎn)化、侵襲和轉(zhuǎn)移等過(guò)程具有重要作用。Warburg等發(fā)現(xiàn)癌細(xì)胞會(huì)比正常細(xì)胞消耗更多的葡萄糖[17]。葡萄糖進(jìn)入細(xì)胞后參與細(xì)胞質(zhì)內(nèi)的糖酵解活動(dòng)，其終產(chǎn)物丙酮酸經(jīng)過(guò)酶促反應(yīng)轉(zhuǎn)化生產(chǎn)乳酸；乳酸可通過(guò)多種機(jī)制促進(jìn)腫瘤的血管形成、細(xì)胞遷移和逃避免疫監(jiān)視[18]。本研究中部分DEGs參與丙酮酸代謝活動(dòng)，這提示癌細(xì)胞可能是通過(guò)糖酵解途徑改變機(jī)體內(nèi)能量代謝方式，從而促進(jìn)COAD的發(fā)生和發(fā)展。此外，本研究中部分DEGs富集于膽汁分泌信號(hào)通路，提示COAD患者的病情進(jìn)展與膽汁分泌具有密切聯(lián)系。已有研究證實(shí)，相對(duì)于健康人群，結(jié)直腸癌患者的腸道微生物組群落發(fā)生變化[19-20]。Ridlon等[21]發(fā)現(xiàn)，飲食中飽和脂肪會(huì)誘導(dǎo)膽汁分泌增多，進(jìn)入腸道后產(chǎn)生脫氧膽酸和石膽酸等二級(jí)膽汁酸，激活細(xì)胞信號(hào)級(jí)聯(lián)反應(yīng)，從而促進(jìn)結(jié)直腸癌細(xì)胞的增殖和遷徙。

鈣激活的氯離子通道蛋白參與細(xì)胞內(nèi)信號(hào)傳導(dǎo)并激活特定的細(xì)胞應(yīng)答，如與癌癥相關(guān)的增殖、凋亡、遷移和血管生成[22]，并被認(rèn)為是新興的藥物靶點(diǎn)[23-24]。CLCA1基因是鈣激活的CLCA家族成員，胃腸道中CLCA1主要在小腸、結(jié)腸和闌尾表達(dá)，為胃腸道提供防護(hù)作用。目前，針對(duì)CLCA1在結(jié)直腸癌中的機(jī)制已有研究。CLCA1可以通過(guò)提高黏附分子E-cadherin和腸堿性磷酸酶的表達(dá)從而促進(jìn)腸上皮分化[25]，而敲除CLCA1(Caco-2細(xì)胞系)則可以抑制細(xì)胞分化并促進(jìn)細(xì)胞增殖[26]。同時(shí)，有研究發(fā)現(xiàn)CLCA1表達(dá)水平的升高可抑制Wnt信號(hào)通路和上皮-間充質(zhì)轉(zhuǎn)化(epithelial-mesenchymal transition，EMT)過(guò)程，表明該基因具有腫瘤抑制作用[25，27-28]。此外，c-myc這一原癌基因的產(chǎn)物參與細(xì)胞增殖和凋亡的調(diào)節(jié)，有研究發(fā)現(xiàn)CLCA1的轉(zhuǎn)錄則與c-myc的轉(zhuǎn)錄具有一定相關(guān)性[29]，但其中具體機(jī)制仍有待明確。

尋找新的治療靶點(diǎn)一直是結(jié)直腸癌研究的關(guān)注要點(diǎn)。目前已有包括西妥昔單抗、帕尼單抗在內(nèi)的以EGFR為靶點(diǎn)的藥物被批準(zhǔn)應(yīng)用于臨床中的靶向治療，但患者仍表現(xiàn)出耐藥性[30]。離子通道的功能表達(dá)及其受類固醇激素和生長(zhǎng)因子的調(diào)控是卵巢癌發(fā)生發(fā)展的重要組成部分，也可能與患者產(chǎn)生耐藥性有關(guān)[31]。Musrap等[32]發(fā)現(xiàn)，CLCA1是聚集形成的卵巢癌細(xì)胞中上調(diào)較多的蛋白；而使用氯離子通道阻滯劑或敲除CLCA1則會(huì)降低癌細(xì)胞形成聚集體的能力，表明CLCA1可能是惡性腫瘤新的治療靶點(diǎn)[33]。CLCA1的低表達(dá)與結(jié)直腸癌的低生存率和高疾病復(fù)發(fā)率有關(guān)[22]，將其作為靶標(biāo)進(jìn)行靶向藥物研發(fā)，有望為COAD的治療提供新方法。

本文的不足之處在于，雖基于WGCNA方法對(duì)COAD潛在的預(yù)后基因展開了一系列的生物信息學(xué)分析和預(yù)測(cè)，但該分析結(jié)果仍需進(jìn)一步的基礎(chǔ)或臨床實(shí)驗(yàn)加以驗(yàn)證。

綜上，進(jìn)一步明確CLCA1對(duì)于COAD發(fā)生和發(fā)展的具體機(jī)制，可能有利于提高臨床預(yù)后的判斷力和個(gè)性化治療方案的優(yōu)化。建議將CLCA1作為靶點(diǎn)進(jìn)行針對(duì)性的臨床研究和靶向藥物研發(fā)，使得更多COAD患者獲益。