溫中慶， 張 強， 孫雨頡， 趙浩楠， 康 寧

天津中醫(yī)藥大學(xué)中西醫(yī)結(jié)合學(xué)院，天津 301617

結(jié)直腸癌(colorectal cancer，CRC)，是發(fā)生于結(jié)直腸部位的消化道惡性腫瘤，結(jié)腸腺癌(colon adenocarcinoma，COAD)是結(jié)直腸癌中最常見的病理類型[1]。據(jù)美國癌癥協(xié)會報道，2019年CRC在男性中發(fā)病率居第2位，在女性中發(fā)病率居第3位。結(jié)直腸癌患者相對5年生存率為65%，而Ⅳ期患者5年生存率僅12%[2]。早期診斷能夠顯著改善結(jié)直腸癌患者的預(yù)后[3]，然而多數(shù)患者發(fā)現(xiàn)時已發(fā)展至中晚期，因此挖掘結(jié)直腸癌發(fā)生發(fā)展過程中的生物標(biāo)志物具有重要意義。

隨著高通量測序技術(shù)的發(fā)展，生物信息學(xué)成為了尋找生物標(biāo)志物的重要手段。加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis，WGCNA)是一種系統(tǒng)生物學(xué)方法，用于檢測高度相關(guān)的基因群組成的基因模塊與特定的臨床性狀之間的相關(guān)性[4]，已經(jīng)被廣泛應(yīng)用于鑒定與多種癌癥臨床特征相關(guān)的基因模塊，包括用于鑒定與結(jié)直腸癌復(fù)發(fā)[5]、微衛(wèi)星不穩(wěn)定[6]等相關(guān)的基因模塊。

GEO(The Gene Expression Omnibus)是一個存儲高通量基因表達(dá)數(shù)據(jù)集的公共數(shù)據(jù)庫[7]。本研究采用NCBI-GEO數(shù)據(jù)庫的GSE89076數(shù)據(jù)集[8]，該數(shù)據(jù)集包含275例結(jié)直腸癌患者的腫瘤組織及其癌旁正常組織的基因表達(dá)譜。TCGA(The Cancer Genome Atlas)數(shù)據(jù)庫是由美國國立衛(wèi)生院啟動的項目，旨在通過對30多種類型人類腫瘤進(jìn)行大規(guī)模基因組測序，提供公開可用的癌癥基因組數(shù)據(jù)，從而推動改進(jìn)診斷標(biāo)準(zhǔn)、治療方法，并最終預(yù)防癌癥[9]。本文采用了來自TCGA-COAD的RNA-seq數(shù)據(jù)，其中包含473個結(jié)直腸癌樣本及41個正常樣本。

本研究基于生物信息學(xué)分析方法研究與結(jié)直腸癌發(fā)生發(fā)展相關(guān)的高風(fēng)險致病基因及其相關(guān)通路，為結(jié)直腸癌的臨床診斷及治療提供潛在靶點。

1 材料與方法

1.1 數(shù)據(jù)下載與預(yù)處理

通過TCGA數(shù)據(jù)庫(https：//portal.gdc.cancer.gov/)下載結(jié)直腸癌樣本和正常樣本的基因表達(dá)數(shù)據(jù)，及其相應(yīng)的臨床數(shù)據(jù)。將Primary site選擇為“colon”，Program選擇為“TCGA”，Project選擇為“TCGA-COAD”，Data category為“transcriptome profiling”，Data type為“Gene expression quantification”，Experimental strategy為“RNA-seq”，Workflow類型為“HTseq-counts”。用R語言處理數(shù)據(jù)，以|logFC|>1，adj.P<0.05為條件篩選差異表達(dá)基因(differentially expressed genes，DEGs)。

通過GEO數(shù)據(jù)庫(https：//www.ncbi.nlm.nih.gov/gds)下載芯片數(shù)據(jù)集GSE89076[8]，芯片平臺為Agilent-039494 SurePrint G3 Human GE v2 8x60K Microarray 039381(Feature Number version)。該微陣列數(shù)據(jù)集包含275例結(jié)直腸癌患者腫瘤組織及其癌旁正常組織的基因表達(dá)數(shù)據(jù)。根據(jù)注釋文件將探針矩陣進(jìn)行注釋，轉(zhuǎn)化為基因矩陣，刪除重復(fù)探針，并以|logFC|>1，adj.P<0.05為條件篩選差異表達(dá)基因(DEGs)。

1.2 加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析

通過R語言中的WGCNA包[4]分別對TCGA-COAD和GSE89076[8]的數(shù)據(jù)構(gòu)建共表達(dá)網(wǎng)絡(luò)。首先，對樣本進(jìn)行聚類，刪除離群樣本；其次，通過篩選軟閾值β，構(gòu)建更符合無尺度網(wǎng)絡(luò)特征的網(wǎng)絡(luò)。通過設(shè)置每個模塊最少基因數(shù)為50，切割高度為0.25進(jìn)行模塊聚類，并對相似模塊進(jìn)行合并，計算模塊顯著性(module significance，MS)。最后計算基因顯著性(gene significance，GS)和模塊身份(module membership，MM)。

1.3 差異共表達(dá)基因GO和KEGG富集及核心基因篩選

用R包VennDiagram[10]將WGCNA篩選出與臨床特征相關(guān)性最強模塊所包含的基因與差異表達(dá)基因取交集，作為差異共表達(dá)基因。隨后通過R包clusterProfiler[11]對差異共表達(dá)基因進(jìn)行GO和KEGG富集分析。隨后采用STRING(https：//www.string-db.org/)在線數(shù)據(jù)庫[12]，選擇score>0.9預(yù)測候選核心基因之間的相互作用，構(gòu)建蛋白質(zhì)互作(protein-protein interaction，PPI)網(wǎng)絡(luò)，用Cytoscape[13]對PPI網(wǎng)絡(luò)進(jìn)行可視化，并通過Cytoscape的CytoHubba插件計算每個節(jié)點的MCC(maximal clique centrality)[14]。在本研究中，MCC值最高的10個基因確定為候選核心基因。為了驗證候選核心基因的預(yù)后價值，我們通過GEPIA2數(shù)據(jù)庫[15]考察了候選核心基因與CRC患者總生存期及無病生存期之間的關(guān)系。以P<0.05為差異具有統(tǒng)計學(xué)意義，將其定義為核心基因，并對核心基因進(jìn)行進(jìn)一步分析。

1.4 人蛋白圖譜數(shù)據(jù)庫驗證生存相關(guān)的核心基因蛋白質(zhì)表達(dá)

人蛋白圖譜數(shù)據(jù)庫(The Human Protein Atlas，HPA，https：//www.proteinatlas.org/)提供了基于免疫組織化學(xué)(immunohistochemistry，IHC)對蛋白質(zhì)相對豐度的檢測[16]。我們通過HPA數(shù)據(jù)庫分析了生存相關(guān)核心基因在結(jié)直腸癌組織和正常組織中的蛋白表達(dá)情況。

1.5 單基因GSEA分析

為了確定核心基因相關(guān)信號通路，我們將TCGA數(shù)據(jù)集中的腫瘤樣本根據(jù)基因表達(dá)中位數(shù)分為高表達(dá)組和低表達(dá)組，選擇c2.cp.kegg.v7.2.symbols.gmt基因集作為參考基因集，使用Gene Set Enrichment Analysis(GSEA)[17](https：//www.gsea-msigdb.org/gsea/index.jsp)軟件v4.1.0進(jìn)行分析，歸一化富集分析得分(NES)是通過分析1000種排列來確定。根據(jù)先前的研究，我們按照|NES|>1、P<0.05和錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)<0.25為標(biāo)準(zhǔn)篩選基因集[18-19]。

1.6 細(xì)胞培養(yǎng)

人正常結(jié)腸上皮細(xì)胞HCoEpiC和結(jié)直腸癌細(xì)胞系HCT116、LoVo均來源于ATCC。HCoEpiC細(xì)胞在含有10%胎牛血清(BI，以色列)和1%青霉素-鏈霉素(Solarbio，中國)的高糖DMEM培養(yǎng)液(biosharp，中國)中培養(yǎng)；HCT116細(xì)胞和LoVo細(xì)胞在含有10%胎牛血清(ExCell，中國)和1%青霉素-鏈霉素(Solarbio，中國)的RPMI-1640培養(yǎng)液(biosharp，中國)中培養(yǎng)。所有細(xì)胞均置于37℃、5% CO2的培養(yǎng)箱中培養(yǎng)。

1.7 Western blot實驗

用RIPA裂解細(xì)胞，12000 r/min離心10 min后吸出上清，通過BCA蛋白定量試劑盒進(jìn)行蛋白定量。將各樣品稀釋至同一濃度，取等量樣品在10%SDS-PAGE中電泳，使用F2RL1(Santa Cruz，1∶100稀釋)及β-actin(中杉金橋，1∶3000稀釋)一抗，置于4℃過夜孵育，加入二抗后置于室溫孵育1 h，使用高敏化學(xué)發(fā)光試劑(Tanon，中國)顯影。使用Image J進(jìn)行光密度分析。

1.8 統(tǒng)計學(xué)方法

采用IBM SPSS 23.0進(jìn)行統(tǒng)計分析，組間首先采用單因素方差分析，隨后通過LSD事后檢驗進(jìn)行多重比較，以P<0.05為差異有統(tǒng)計學(xué)意義。采用GraphPad Prism version 7.0進(jìn)行統(tǒng)計作圖。

2 結(jié)果

2.1 TCGA-COAD及GSE89076差異基因表達(dá)譜構(gòu)建

通過R語言limma包[20]對TCGA數(shù)據(jù)庫下載結(jié)直腸癌樣本和正常樣本及GSE89076數(shù)據(jù)集進(jìn)行分析，并以|logFC|>1，P<0.05為條件篩選DEGs，并通過使用R語言ggplot2[21]繪制火山圖和熱圖，其中紅色代表表達(dá)上調(diào)的基因，藍(lán)色代表表達(dá)下調(diào)的基因。在TCGA-COAD中共得到3481個DEGs，其中包括1208個上調(diào)基因和2273個下調(diào)基因(圖1)。在GSE89076數(shù)據(jù)集中共得到3460個DEGs，其中包括1753個上調(diào)基因和1707個下調(diào)基因(圖2)。

A：DEGs火山圖；B：DEGs熱圖圖1 TCGA-COAD差異表達(dá)基因(DEGs)分析Fig.1 Analysis of differentialy expressed genes(DEGs)in TCGA-COAD

A：DEGs火山圖；B：DEGs熱圖圖2 GSE89076差異表達(dá)基因(DEGs)分析Fig.2 Analysis of differentialy expressed genes(DEGs)in GSE89076

2.2 TCGA-COAD及GSE89076加權(quán)基因共表達(dá)模塊構(gòu)建

通過R語言的WGCNA包[4]構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)。對樣本進(jìn)行聚類，去除離群值后繪制樣本聚類樹。在TCGA-COAD中選擇β=3(R2=0.9)構(gòu)建無標(biāo)度網(wǎng)絡(luò)(圖3)；在GSE89076數(shù)據(jù)集中選擇β=11(R2=0.9)構(gòu)建無標(biāo)度網(wǎng)絡(luò)(圖4)。本研究共確定了TCGA-COAD中的11個共表達(dá)模塊及GSE89076數(shù)據(jù)集的22個共表達(dá)模塊(其中包含1個無法被歸類的基因組成的grey模塊)。接下來繪制模塊-臨床特征相關(guān)性熱圖，以評估模塊與臨床特征(腫瘤與正常)之間的相關(guān)性。其中TCGA-COAD中yellow模塊與腫瘤組織相關(guān)性最強(r=-0.87，P<0.05)；GSE89076數(shù)據(jù)集中blue模塊與腫瘤組織相關(guān)性最強(r=-0.77，P<0.05)。因此選擇這2個模塊為重要模塊進(jìn)行進(jìn)一步分析。

A：樣本聚類及臨床特征熱圖；B：通過無標(biāo)度擬合指數(shù)和平均連通性篩選軟閾值；C：共表達(dá)網(wǎng)絡(luò)模塊聚類樹狀圖；D：基因模塊與臨床特征相關(guān)性熱圖圖3 確定TCGA-COAD中與臨床特征相關(guān)的模塊Fig.3 Identification of modules associated with clinical features in the TCGA-COAD

A：樣本聚類及臨床性狀熱圖；B：通過無標(biāo)度擬合指數(shù)和平均連通性篩選軟閾值；C：基因聚類樹狀圖，不同顏色代表不同基因模塊；D：基因模塊與臨床特征相關(guān)性熱圖圖4 確定GSE89076數(shù)據(jù)集中與臨床特征相關(guān)的模塊Fig.4 Identification of modules associated with clinical features in the GSE89076

2.3 差異共表達(dá)基因鑒定及其富集分析

將TCGA-COAD中yellow模塊，GSE89076數(shù)據(jù)集中blue模塊，TCGA-COAD的DEGs和GSE89076數(shù)據(jù)集的DEGs取交集，識別到436個重疊基因作為差異共表達(dá)基因(圖5A)。為了進(jìn)一步探究436個差異共表達(dá)基因的潛在功能，我們使用clusterProfiler包[11]對其進(jìn)行了富集分析。KEGG富集分析主要富集在補體系統(tǒng)、細(xì)胞因子-細(xì)胞因子受體的相互作用途徑、腫瘤壞死因子信號通路、細(xì)胞粘附分子途徑和谷胱甘肽代謝途徑等信號通路(圖5B)。通過GO富集分析，其生物學(xué)過程(biological process，BP)主要在脂質(zhì)分解和脂肪酸代謝過程中富集；細(xì)胞成分(cellular component，CC)分析主要涉及細(xì)胞頂端和頂端質(zhì)膜；分子功能(molecular function，MF)分析中主要富集在陰離子跨膜轉(zhuǎn)運體活性和磷酸酯水解酶活性(圖5C)。

A：DEGs與重要共表達(dá)模塊之間的韋恩圖，共436個交集基因作為差異共表達(dá)基因；B：差異共表達(dá)基因的KEGG富集分析；C：差異共表達(dá)基因的GO富集分析圖5 確定差異共表達(dá)基因及其GO、KEGG富集分析Fig.5 Identification of differentially co-expressed genes and their results from GO and KEGG analysis

2.4 候選核心基因篩選

隨后采用STRING數(shù)據(jù)庫[12]構(gòu)建436個基因的蛋白質(zhì)-蛋白質(zhì)互作(PPI)網(wǎng)絡(luò)，通過Cytoscape[13]對網(wǎng)絡(luò)進(jìn)行可視化，并采用CytoHubba插件的MCC算法從PPI網(wǎng)絡(luò)中篩選排名前10的基因作為候選核心基因[14]：BDKRB2、GCG、EDN2、EDN3、P2RY2、P2RY1、F2RL1、GNA11、SST、ADRA2(圖6)。

紅色節(jié)點代表MCC得分較高的基因，黃色節(jié)點代表MCC得分較低的基因圖6 PPI網(wǎng)絡(luò)及候選核心基因可視化Fig.6 Visualization of the PPI network and the candidate hub genes

2.5 核心基因篩選及其mRNA和蛋白質(zhì)表達(dá)模式驗證

我們用GEPIA2數(shù)據(jù)庫[15]考察了10個候選核心基因的預(yù)后價值。Kaplan-Meier分析顯示僅F2RL1的低表達(dá)與較差的總生存期(overall survival，OS)相關(guān)(圖7)，而低表達(dá)F2RL1與無病生存期(disease-free survival，DFS)無顯著相關(guān)性(圖8)。隨后，通過GEPIA2數(shù)據(jù)庫[15]，對F2RL1在人體正常組織樣本和COAD樣本中的mRNA表達(dá)差異進(jìn)行分析，發(fā)現(xiàn)F2RL1的mRNA水平在腫瘤組織中顯著降低(圖9)。此外，分析HPA數(shù)據(jù)庫[16]，顯示CRC腫瘤組織中F2RL1的蛋白質(zhì)水平低于正常組織(圖10)。上述研究證實，F(xiàn)2RL1在CRC組織中表達(dá)下調(diào)，且F2RL1的低表達(dá)與CRC患者總生存期較短相關(guān)。

A：ADRA2；B：BDKRB2；C：EDN2；D：EDN3；E：F2RL1；F：GCG；G：GNA11；H：P2RY1；I：P2RY2；J：SST圖7 通過GEPIA2在線工具分析10個候選核心基因與COAD患者總生存期的關(guān)系Fig.7 Analysis of relationship between 10 candidate hub genes and OS of COAD patients by using the GEPIA2 online tool

A：ADRA2；B：BDKRB2；C：EDN2；D：EDN3；E：F2RL1；F：GCG；G：GNA11；H：P2RY1；I：P2RY2；J：SST圖8 通過GEPIA2在線工具分析10個候選基因與COAD患者無病生存期的關(guān)系Fig.8 Analysis of relationship between 10 candidate hub genes and DFS of COAD patients by using the GEPIA2 online tool

*P<0.05圖9 使用GEPIA2對TCGA數(shù)據(jù)庫的275個COAD樣本和41個正常樣本中F2RL1 mRNA水平進(jìn)行差異表達(dá)分析Fig.9 Analysis of differential expression of F2RL1 mRNA in 275 COAD samples and 41 normal samples in the TCGA database by using GEPIA2

A:F2RL1在腫瘤樣本中的表達(dá)；B:F2RL1在正常樣本中的表達(dá)圖10 HPA數(shù)據(jù)庫中免疫組織化學(xué)染色顯示F2RL1蛋白在CRC組織和正常組織中的表達(dá)(×100)Fig.10 Immunohistochemistry analysis of F2RL1 in CRC tissues and normal tissues from HPA database(×100)

2.6 對F2RL1進(jìn)行單基因GSEA分析

為了進(jìn)一步確定核心基因的潛在功能，本研究對TCGA-COAD來源的數(shù)據(jù)進(jìn)行了基于KEGG基因集的單基因GSEA。單基因GSEA揭示了F2RL1的上調(diào)與核苷酸切除修復(fù)信號通路(NES=1.77，P<0.01)、P53信號通路(NES=1.70，P=0.009)、蛋白酶體信號通路(NES=1.60，P=0.034)、錯配修復(fù)信號通路(NES=1.57，P=0.037)等呈正相關(guān)(圖11)。

圖11 基于TCGA-COAD數(shù)據(jù)對F2RL1進(jìn)行單基因GSEA分析(僅列出5個最具有代表性的常見KEGG通路)Fig.11 Single-gene GSEA analysis of F2RL1 based on TCGA-COAD data(only top five representative KEGG pathways are listed)

2.7 F2RL1在人結(jié)直腸細(xì)胞系中的表達(dá)

通過Western blot驗證了F2RL1在人正常結(jié)腸上皮細(xì)胞HCoEpiC和結(jié)直腸癌細(xì)胞HCT116、LoVo中的表達(dá)水平。與生信分析預(yù)測結(jié)果一致，結(jié)直腸癌細(xì)胞中F2RL1的蛋白表達(dá)較正常結(jié)腸上皮細(xì)胞中的蛋白表達(dá)顯著下調(diào)(圖12)。

與HCoEpiC比較，*P<0.05圖12 F2RL1蛋白質(zhì)表達(dá)水平驗證Fig.12 Verification of expression levels of F2RL1 protein

3 討論

CRC是最常見的胃腸道腫瘤類型[22]，在中國的發(fā)病率呈上升趨勢[23]。超過90%的CRC由結(jié)腸息肉發(fā)展而來[24]，如果在息肉階段治愈，防止癌變，5年生存率可達(dá)90%以上[25]，但是晚期患者的5年生存率<10%。目前，結(jié)腸鏡和組織活檢是主要用于早期發(fā)現(xiàn)CRC的手段，然而結(jié)腸鏡是一種侵入性手術(shù)，成本高，且需要較高水平的專業(yè)技術(shù)。因此，需要尋找更好的CRC特異性預(yù)后和進(jìn)展生物標(biāo)志物及治療靶點。在本研究中，我們使用WGCNA和差異表達(dá)基因分析，在TCGA-COAD和GSE89076數(shù)據(jù)集中鑒定了436個表達(dá)趨勢相同的差異共表達(dá)基因。進(jìn)一步通過PPI網(wǎng)絡(luò)確定了BDKRB2、GCG、EDN2、EDN3、P2RY2、P2RY1、F2RL1、GNA11、SST、ADRA2這10個候選核心基因。將10個候選基因進(jìn)行預(yù)后價值驗證，最終確定F2RL1為核心基因。

F2RL1(F2R like trypsin receptor 1)編碼的蛋白為蛋白酶激活受體2(proteinase-activated receptor 2，PAR2)，屬于G蛋白偶聯(lián)受體家族[26]。在以往的研究中發(fā)現(xiàn)F2RL1與炎癥[27]、疼痛[28]、膽固醇與脂質(zhì)代謝[29]、肥胖[30]、血管生成等相關(guān)。除此之外，據(jù)文獻(xiàn)報道，F(xiàn)2RL1的mRNA和蛋白表達(dá)在炎癥性腸病和結(jié)直腸癌樣本中均下調(diào)，且F2RL1下調(diào)與較短的生存期相關(guān)[31]。有趣的是，Kawaguchi等[32]發(fā)現(xiàn)F2RL1促進(jìn)腸道腫瘤的發(fā)生；也有研究表明在肺腺癌中低表達(dá)F2RL1與較好的預(yù)后相關(guān)[33]；在肝細(xì)胞癌中F2RL1促進(jìn)遷移和侵襲[34]。本研究通過對F2RL1進(jìn)行單基因GSEA分析，確定了核苷酸切除修復(fù)信號通路、P53信號通路、蛋白酶體信號通路、錯配修復(fù)信號通路在高表達(dá)F2RL1組中被上調(diào)。核苷酸切除修復(fù)(NER)是保護(hù)細(xì)胞免受基因組損傷的主要DNA修復(fù)途徑之一，這條通路的中斷會促進(jìn)癌癥的發(fā)展[35]，并且癌細(xì)胞在參與DNA修復(fù)的相關(guān)基因中顯示出各種類型的突變和異常表達(dá)。這些改變會導(dǎo)致基因組不穩(wěn)定并促進(jìn)癌癥發(fā)生和癌癥進(jìn)展過程。因此DNA修復(fù)中的這些缺陷也被認(rèn)為是癌癥治療的合適靶點[36]；轉(zhuǎn)錄因子P53主要通過調(diào)控其下游凋亡、細(xì)胞周期、DNA修復(fù)、衰老等相關(guān)靶基因轉(zhuǎn)錄，在抑制腫瘤中發(fā)揮極其重要的作用，其突變或失活是幾乎所有腫瘤的標(biāo)志[37-38]；泛素蛋白酶體通過調(diào)控細(xì)胞中大多數(shù)蛋白質(zhì)的降解，從而參與許多細(xì)胞過程的調(diào)控[39]，其功能障礙與多種腫瘤的發(fā)生相關(guān)[40]；錯配修復(fù)系統(tǒng)是維持DNA穩(wěn)定性的保證[41]，DNA錯配修復(fù)缺陷會導(dǎo)致高突變表型，從而導(dǎo)致腫瘤的發(fā)生[42]。這些結(jié)果進(jìn)一步證實了F2RL1可能在CRC發(fā)生發(fā)展中通過上述通路發(fā)揮重要作用。

綜上所述，通過WGCNA和差異表達(dá)基因分析，我們的研究確定了CRC中具有預(yù)后價值的核心基因F2RL1，并且進(jìn)一步驗證了F2RL1在人結(jié)直腸癌細(xì)胞系中表達(dá)下調(diào)。雖然我們通過對多個數(shù)據(jù)庫的CRC相關(guān)數(shù)據(jù)進(jìn)行了生物信息學(xué)分析，但識別核心基因只是篩選CRC生物標(biāo)志物及治療靶點的第一步，我們?nèi)孕枰ㄟ^一系列體外體內(nèi)的實驗驗證F2RL1在CRC發(fā)生發(fā)展和預(yù)后中的價值，以及其發(fā)揮作用的分子機制。