基于加權基因共表達網絡分析篩選結腸腺癌預后關鍵基因

李卓陽， 張皓旻， 劉格良， 陳浩然， 智 鵬， 陳熙勐， 盧學春， 賀培鳳

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

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

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

本研究應用WGCNA方法，對來自TCGA和GEO數據庫的COAD RNA-seq數據，構建COAD的差異基因共表達網絡，篩選與COAD發生發展密切相關的基因模塊，隨后結合Kaplan-Meier方法鑒定COAD預后基因，為COAD預后相關基因的進一步基礎及臨床研究提供方向和指導。

1 材料與方法

1.1 數據獲取從TCGA數據庫[10]中獲取COAD轉錄組數據和相應的臨床信息。其中，轉錄組數據包括398例COAD樣本和39例正常結腸組織樣本；臨床信息包括COAD患者的生存狀態和生存時間。

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

1.2 差異表達基因分析采用R語言軟件包limma，分別對來自TCGA和GEO的轉錄組數據進行數據標準化和差異表達基因(differentially expressed genes，DEGs)分析，得到兩組DEGs。差異表達基因的篩選標準，TCGA為|logFC|≥1、FDR<0.05；GEO為|logFC|≥1、FDR<0.05。當多個探針與一個相同的基因匹配時，則以平均值作為該基因的表達值。最后采用R語言軟件包ggplot2對DEGs的表達模式進行可視化。

1.3 加權基因共表達網絡分析采用R語言軟件包WGCNA，分別對來自TCGA和GEO的轉錄組數據進行加權基因共表達網絡分析。填補缺失值后，構建鄰接矩陣并轉換為拓撲重疊矩陣。采用動態剪切法構建基因聚類樹狀圖，將表達相似的基因聚類為不同的基因共表達模塊，規定每個模塊最少基因數目為50。隨后采用Pearson相關分析計算每個基因與各個模塊、不同臨床特征基因顯著性的相關系數r值和P值，獲取腫瘤相關性最高的模塊及該模塊中的基因。最后采用R語言軟件包VennDiagram，將相關系數最高模塊中的基因與TCGA、GEO的DEGs取交集，得到最終的DEGs，用以后續分析。

1.4 GO富集和KEGG通路富集分析采用R語言軟件包clusterProfiler，對上述DEGs進行GO[13-14]富集分析和KEGG通路[15]富集分析。以FDR<0.05、P<0.05為閾值，篩選富集的GO條目和KEGG通路。

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

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

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

2 結果

2.1 差異表達基因篩選從TCGA數據庫下載COAD數據集，經分析得到差異表達基因3 544個，其中上調基因1 293個，下調基因2 251個。從GEO數據庫下載GSE110224數據集，經分析得到差異表達基因515個，其中上調基因223個，下調基因292個(見圖1)。

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

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

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

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

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

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

注：橫坐標代表富集在各GO條目的基因數目比率，縱坐標代表GO條目名稱。

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

注：橫坐標代表富集在各通路的基因數目比率，縱坐標代表通路名稱。

2.4 PPI網絡構建和COAD預后相關核心基因篩選PPI網絡中共出現153個節點和385條連線(見圖7A)。根據Closeness算法，凝聚素Ⅰ復合物亞基G(non-SMC condensin I complex subunit G，NCAPG)、細胞分裂周期6(cell division cycle 6 homolog，CDC6)、核受體亞家族1，組H，成員4(nuclear receptor subfamily 1, group H, member 4，NR1H4)、甲狀腺素受體結合因子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)、蔗糖酶異麥芽糖酶復合物(sucrase isomaltase，SI)、周期蛋白B1(Cyclin B1，CCNB1)、叉頭框M1(forkhead box M1，FOXM1)等10個基因為其中的核心基因(見圖7B)。

注：A：153個DEGs的PPI網絡；B：153個DEGs的核心基因。

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

圖8 OS分析結果

圖9 10個核心基因的DFS分析結果

注：A;正常結腸組織，B;結腸腺癌組織。

3 討論

本研究通過對TCGA和GEO數據庫中的COAD轉錄組數據進行差異表達基因分析和加權基因共表達網絡分析，發現COAD主要與患者體內的153個基因異常表達有關。富集分析發現這些基因大多與離子運輸、激素代謝等生命活動有關。最后通過蛋白互作網絡分析和生存分析，發現CLCA1基因與COAD患者的不良預后顯著相關。

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

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

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

本文的不足之處在于，雖基于WGCNA方法對COAD潛在的預后基因展開了一系列的生物信息學分析和預測，但該分析結果仍需進一步的基礎或臨床實驗加以驗證。

綜上，進一步明確CLCA1對于COAD發生和發展的具體機制，可能有利于提高臨床預后的判斷力和個性化治療方案的優化。建議將CLCA1作為靶點進行針對性的臨床研究和靶向藥物研發，使得更多COAD患者獲益。