共表達(dá)網(wǎng)絡(luò)鑒定食管鱗癌與腺癌關(guān)鍵差異表達(dá)分子

李孟祥, 劉軻, 阮豪杰, 高社干, 齊義軍*, 馮笑山

(1.河南科技大學(xué) 信息工程學(xué)院，河南 洛陽(yáng) 471023；2.河南科技大學(xué) 臨床醫(yī)學(xué)院/第一附屬醫(yī)院/河南省腫瘤表觀遺傳重點(diǎn)實(shí)驗(yàn)室，河南 洛陽(yáng) 471003)

食管癌(esophageal cancer, EC)位居中國(guó)惡性腫瘤發(fā)病率和死亡率的第4位，中國(guó)每年確診的EC病例和死亡病例約占全球總確診病例和死亡病例的一半[1].2015年，中國(guó)EC的新發(fā)病例、新增死亡病例約為24.6萬(wàn)例、18.8萬(wàn)例[2].EC主要有兩種組織病理學(xué)類型：食管鱗狀細(xì)胞癌(esophageal squamous cell carcinoma, ESCC)和食管腺癌(esophageal adenocarcinoma, EAC).在中國(guó)，至少90%的EC為ESCC，而歐美等西方國(guó)家卻以EAC為主[3-4].

ESCC大多發(fā)生在食管中、上段，病因?qū)W上與抽煙、酗酒等密切相關(guān).EAC主要發(fā)生于食管下段與賁門交界處，病因?qū)W上與肥胖、胃食管反流等關(guān)系密切.此外，ESCC與肺、宮頸等組織部位的鱗癌基因突變模式類似[5]，而EAC和胃腺癌分子病變模式相近.因此，鑒定起源于食管但病理類型不同的ESCC與EAC特異性分子表達(dá)譜，是靶向治療和個(gè)體化治療的基礎(chǔ).

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighted gene co-expression network analysis, WGCNA)通過(guò)構(gòu)建共表達(dá)網(wǎng)絡(luò)確定基因模塊，分析模塊與性狀表型之間的相關(guān)性，從而挖掘與臨床表型相關(guān)的關(guān)鍵分子[6].本研究應(yīng)用公共數(shù)據(jù)庫(kù)中EC的mRNA表達(dá)譜測(cè)序數(shù)據(jù)，鑒定EAC和ESCC差異表達(dá)基因，然后應(yīng)用WGCNA算法構(gòu)建基因共表達(dá)網(wǎng)絡(luò)模塊，確定模塊內(nèi)關(guān)鍵基因以鑒別EAC和ESCC并進(jìn)行驗(yàn)證．

1 材料和方法

1.1 EC的mRNA表達(dá)譜數(shù)據(jù)

從癌癥基因圖譜數(shù)據(jù)庫(kù)(the cancer genome atlas, TCGA)和基因表達(dá)數(shù)據(jù)庫(kù)(gene expression omnibus, GEO)下載EC的mRNA轉(zhuǎn)錄組數(shù)據(jù)及相應(yīng)臨床信息，分別為TCGA-EC和GSE26866.將表達(dá)譜數(shù)據(jù)與臨床信息匹配后，共納入來(lái)源于TCGA-EC的78例EAC和80例ESCC；GSE26866數(shù)據(jù)集包括20例Barrett’s食管(Barrett’s esophagus, BE)、21例EAC、9例ESCC和19例食管鱗狀上皮(esophageal squamous epithelium, ESE)樣本(表1)．

表1 入組數(shù)據(jù)集的組織類型和樣本數(shù)量

1.2 TCGA-EC差異表達(dá)基因鑒定

用R語(yǔ)言中DESeq2程序包處理TCGA-EC mRNA表達(dá)譜數(shù)據(jù)，差異表達(dá)基因篩選條件為：差異倍數(shù)(fold change, FC)>2，錯(cuò)誤發(fā)現(xiàn)率(false discovery rate, FDR)<0.05且base mean >10；然后將所有基因表達(dá)值進(jìn)行歸一化處理，確定適合樣本相關(guān)性計(jì)算的基因表達(dá)譜數(shù)據(jù).

1.3 WGCNA網(wǎng)絡(luò)模塊構(gòu)建

在WGCNA構(gòu)建的基因共表達(dá)網(wǎng)絡(luò)中，同一個(gè)模塊內(nèi)的基因表達(dá)模式相似，而不同模塊的基因表達(dá)模式差別較大，通過(guò)加權(quán)使得整個(gè)網(wǎng)絡(luò)的連接服從無(wú)標(biāo)度網(wǎng)絡(luò)分布[7].本研究在構(gòu)建無(wú)標(biāo)度網(wǎng)絡(luò)時(shí)，確定的軟閾值β須滿足：R2>0.9，且平均連接度小于100.通過(guò)名為“WGCNA”的程序包構(gòu)建基因聚類樹，并采用動(dòng)態(tài)剪切算法鑒定基因模塊[8].

根據(jù)模塊第一主成分Module Eigengenes(ME)計(jì)算模塊內(nèi)所有基因與ME的相關(guān)性，確定基因的模塊隸屬度(module membership, MM)；根據(jù)ME與EAC、ESCC表型的相關(guān)性，確定基因顯著性(gene significance, GS)；MM與GS均滿足一定條件的基因?yàn)槟K的節(jié)點(diǎn)基因．

1.4 富集分析及獨(dú)立驗(yàn)證

對(duì)重點(diǎn)模塊的節(jié)點(diǎn)基因進(jìn)行GO和KEGG富集分析，其中GO富集分析包括細(xì)胞成分(cellular component, CC)、生物學(xué)過(guò)程(biological process, BP)和分子功能(molecular function, MF)；KEGG富集分析則是在所有已知的生物學(xué)通路(pathway)中進(jìn)行.根據(jù)超幾何分布檢驗(yàn)顯著性FDR，確定模塊內(nèi)顯著富集的生物學(xué)功能和通路.使用ClusterProfiler程序包進(jìn)行富集分析及相關(guān)繪圖[9].

應(yīng)用3組獨(dú)立數(shù)據(jù)驗(yàn)證上述關(guān)鍵基因鑒別ESCC和EAC的準(zhǔn)確性，分別是GSE26866、肺鱗狀細(xì)胞癌(lung squamous cell carcinoma, LUSC)和肺腺癌(lung adenocarcinoma, LUAD)、癌細(xì)胞數(shù)據(jù)庫(kù)(cancer cell line encyclopedia, CCLE)中23個(gè)ESCC細(xì)胞系和1個(gè)EAC細(xì)胞系.

1.5 統(tǒng)計(jì)學(xué)分析

用R語(yǔ)言(3.6.3版本)進(jìn)行統(tǒng)計(jì)學(xué)分析.兩組樣本表達(dá)譜差異比較，采用Wilcoxon檢驗(yàn)；多組樣本表達(dá)譜差異比較，采用Kruskal-Wallis檢驗(yàn)，并進(jìn)行P值的FDR校正.以P<0.05或FDR<0.05為差異有統(tǒng)計(jì)學(xué)意義.

2 結(jié)果

2.1 EAC和ESCC差異表達(dá)基因

根據(jù)本研究確定的差異表達(dá)基因篩選標(biāo)準(zhǔn)，分析TCGA數(shù)據(jù)集中EAC和ESCC基因表達(dá)數(shù)據(jù)，共鑒定到EAC和ESCC的差異表達(dá)基因3745個(gè)，ESCC中上調(diào)和下調(diào)的基因個(gè)數(shù)分別為1 885個(gè)(50.33%)和1 860個(gè)(49.67%，圖1).

紅色和綠色分別代表食管鱗癌中上調(diào)和下調(diào)的基因Red and green represent up-regulated and down-regulated genes in esophageal squamous cell carcinoma, respectively圖1 TCGA數(shù)據(jù)集中食管鱗癌與食管腺癌差異表達(dá)基因火山圖Fig.1 Volcano map of DEGs in esophageal squamous cell carcinoma and esophageal adenocarcinoma in the TCGA dataset

2.2 共表達(dá)網(wǎng)絡(luò)模塊分析

基于ESCC和EAC的3745個(gè)差異表達(dá)基因，對(duì)軟閾值β取1～20的網(wǎng)絡(luò)拓?fù)溥M(jìn)行分析，在標(biāo)度獨(dú)立性和平均連通性條件下選擇β=5，此時(shí)網(wǎng)絡(luò)無(wú)標(biāo)度拓?fù)鋽M合指數(shù)大于0.9且平均連接度小于100，網(wǎng)絡(luò)結(jié)構(gòu)近似于無(wú)標(biāo)度網(wǎng)絡(luò).通過(guò)動(dòng)態(tài)剪切樹算法將相似性大于0.8的模塊合并為同一個(gè)模塊，最終共得到18個(gè)模塊，其余不能歸類到上述18個(gè)模塊的304個(gè)基因歸為灰色(grey)模塊(圖2).

上方為基因聚類圖譜；下方為動(dòng)態(tài)剪切后基因所屬模塊，不同顏色代表不同模塊.The upper part is the gene clustering map and the bottom is the module to which the gene belongs after dynamic shearing with different colors denoting different modules圖2 食管鱗癌和食管腺癌差異表達(dá)基因和共表達(dá)網(wǎng)絡(luò)模塊分布Fig.2 Distribution of differentially expressed genes between esophageal squamous cell carcinoma and esophageal adenocarcinoma, and co-expression network modules by WGCNA

應(yīng)用相關(guān)性分析計(jì)算19個(gè)模塊中每個(gè)模塊ME與ESCC/EAC、BMI、性別、年齡和臨床分期的相關(guān)性，圖3顯示19個(gè)模塊在4種不同表型中的分布.與其他3個(gè)表型特征相比，不同組織學(xué)類型食管腫瘤(ESCC vs EAC)與模塊的相關(guān)性最強(qiáng)，其中9個(gè)模塊與ESCC呈正相關(guān)，10個(gè)模塊與EAC呈正相關(guān).19個(gè)模塊與BMI的相關(guān)性僅次于組織學(xué)類型，令人感興趣的是，19個(gè)模塊在食管腫瘤的組織學(xué)類型與BMI中的分布完全相反，與ESCC呈正相關(guān)的模塊則與BMI呈負(fù)相關(guān)，反之亦然.19個(gè)模塊在年齡特征和臨床分期特征中分布與BMI相似，但與ESCC/EAC呈負(fù)相關(guān)，但相關(guān)性低于BMI.所有模塊與性別相關(guān)性系數(shù)均介于-0.2到0.11之間，相關(guān)性較小(圖3).

與EAC呈正相關(guān)的10個(gè)模塊中，Blue模塊的相關(guān)性最大(cor =-0.95，P=7e-79)，該模塊共包含1 265個(gè)基因，這些基因在EAC中高表達(dá)；而Turquoise模塊與ESCC呈正相關(guān)(cor=0.89，P=9e-55)，該模塊中的1 454個(gè)基因在ESCC中高表達(dá).4a和4b分別顯示了Blue和Turquoise模塊內(nèi)所有基因的模塊隸屬度MM與基因顯著性GS的分布關(guān)系.

圖3 模塊與ESCC/EAC、BMI、Gender、Age、 Stage(腫瘤臨床分期)的相關(guān)性熱圖Fig.3 Correlations between heat map of modules and ESCC/EAC, BMI, Gender, Age, Stage

2.3 EAC和ESCC代表性基因確定及驗(yàn)證

分別計(jì)算Blue模塊和Turquoise模塊中構(gòu)成基因的MM與GS之和，以確定每個(gè)模塊中的節(jié)點(diǎn)基因.結(jié)果表明，EPS8L3和SOX15分別是Blue模塊和Turquoise模塊中MM與GS之和最大的基因，即最重要的節(jié)點(diǎn)基因，因此確定EPS8L3和SOX15分別為EAC和ESCC的代表性基因.

為驗(yàn)證本實(shí)驗(yàn)確定的EAC和ESCC代表性基因，本實(shí)驗(yàn)應(yīng)用3組不同來(lái)源的樣本進(jìn)行驗(yàn)證.在TCGA的80例ESCC和78例EAC樣本中，EPS8L3在EAC和ESCC中的表達(dá)值分別為12.83±1.21、5.87±0.86(P<0.000 1，圖4c)，SOX15的表達(dá)值分別為7.91±1.47、12.34±1.16(P<0.000 1，圖4d).在GSE26866數(shù)據(jù)集中，EPS8L3在21例EAC中的表達(dá)明顯高于ESCC(校正后P<0.001，圖4e)，而9例ESCC組織中SOX15的表達(dá)明顯高于EAC(校正后P<0.001，圖4f)，EPS8L3和SOX15在ESCC中的表達(dá)值分別是-3.98±0.82、3.15±1.99.同樣，EPS8L3在EAC細(xì)胞系OE19中表達(dá)值為48.89，顯著高于23株ESCC細(xì)胞系中的平均表達(dá)值0.059；相反，SOX15在23株ESCC細(xì)胞系中平均表達(dá)值明顯高于OE19細(xì)胞系(25.72和0.335).為進(jìn)一步驗(yàn)證EPS8L3和SOX15是腺上皮和鱗狀上皮的特異性標(biāo)志物分子，本實(shí)驗(yàn)應(yīng)用TCGA數(shù)據(jù)庫(kù)中的498例LUSC和476例LUAD做進(jìn)一步驗(yàn)證，EPS8L3和SOX15在LUSC和LUAD中的表達(dá)模式與ESCC和EAC一致(P<0.01，圖4g；P<0.01，圖4h).上述實(shí)驗(yàn)結(jié)果表明，EPS8L3和SOX15分別是上皮來(lái)源的腺癌和鱗狀細(xì)胞癌特異性分子，可用于二者的鑒別診斷.

(a-b): blue和turquoise模塊內(nèi)基因的MM與GS特征散點(diǎn)圖； (c-d): EPS8L3和SOX15在數(shù)據(jù)集TCGA中的表達(dá)值比較； (e-f): EPS8L3和SOX15在數(shù)據(jù)集GSE26886中的表達(dá)譜比較，差異顯著性為FDR法矯正后的P值； (g-h): EPS8L3和SOX15在數(shù)據(jù)集LUSC、LUAD中的表達(dá)譜比較. 兩組樣本比較用Wilcoxon檢驗(yàn)；多組樣本比較用Kruskal-Wallis檢驗(yàn)；1)P<0.001.(a-b): Scatter plots of MM and GS characteristics of genes in blue and turquoise module; (c-d): Comparison of expression profiles of EPS8L3 and SOX15 in the data sets ESCC and EAC; (e-f): Comparison of expression profiles of EPS8L3 and SOX15 in data set GSE26886, the P value was corrected by the FDR method; (g-h): Comparison of expression profiles of EPS8L3 and SOX15 in data sets LUSC and LUAD. The Wilcoxon test was used for the comparison of two groups of samples; the Kruskal-Wallis test was used for the comparison of multiple groups of samples; 1) P<0.001.圖4 關(guān)鍵節(jié)點(diǎn)基因表達(dá)值比較Fig.4 Comparation of key hub genes from modules and their expression values.

2.4 基因生物學(xué)功能和通路富集分析

以MM>0.8和GS>0.7為標(biāo)準(zhǔn)篩選模塊節(jié)點(diǎn)基因，Blue模塊和Turquoise模塊分別鑒定了259個(gè)和66個(gè)節(jié)點(diǎn)基因.對(duì)Blue模塊中259個(gè)基因做GO富集分析，結(jié)果表明這些基因的功能主要富集于細(xì)胞級(jí)性(Apical part of cell)、肌動(dòng)蛋白相關(guān)的細(xì)胞形態(tài)(Actin-based cell projection)等(圖5a).對(duì)Turquoise模塊66個(gè)節(jié)點(diǎn)基因的GO和KEGG富集分析顯示，這些基因主要富集于表皮發(fā)育(epidermis development)、角質(zhì)細(xì)胞分化(keratinocyte differentiation)、表皮細(xì)胞分化(epidermal cell differentiation)、雌激素信號(hào)通路(estrogen signaling pathway)、Rap1信號(hào)通路(rap1 signaling pathway)等生物學(xué)功能(圖5b).

(a): Blue模塊GO富集分析結(jié)果； (b): Turquoise模塊的GO富集分析結(jié)果.(a): GO enrichment analysis result of blue module; (b): GO enrichment analysis result of turquoise module.圖5 Blue和Turquoise模塊組成基因富集分析Fig.5 Enrichment analysis of blue and turquoise

3 討論

在病理類型、病因?qū)W、好發(fā)部位等方面，中國(guó)與歐美等西方國(guó)家EC完全不同[10]，因此，我國(guó)以ESCC為主的臨床治療方案不能完全借鑒歐美等國(guó)家EAC為主的治療方案，亟需制定適合中國(guó)EC，尤其是ESCC的治療原則及有效治療方案.

應(yīng)用芯片及高能量測(cè)序技術(shù)和生物信息學(xué)方法鑒定不同腫瘤、同一部位不同組織類型的腫瘤，能夠?yàn)槟[瘤的精準(zhǔn)治療提供有效的分子靶點(diǎn).有研究表明[11]，SOX2、PIK3CA、MYC、CCND1、FGR1、GATA4、GATA6等基因位點(diǎn)擴(kuò)增在ESCC和EAC中明顯不同，含有SOX2和PIK3CA基因位點(diǎn)的3q在ESCC中擴(kuò)增頻率顯著高于EAC中(60% vs 15%).應(yīng)用外顯子測(cè)序數(shù)據(jù)分析在美國(guó)ESCC和EAC中的熱點(diǎn)基因的突變頻率，發(fā)現(xiàn)NOTCH1失活突變?cè)贓SCC中存在較高的突變頻率(21%)，而EAC中不存在NOTCH1的失活突變.進(jìn)一步比較美國(guó)和中國(guó)ESCC熱點(diǎn)基因分子差異，發(fā)現(xiàn)美國(guó)ESCC中NOTCH1的失活突變頻率遠(yuǎn)遠(yuǎn)高于中國(guó)的ESCC(21% vs 2%)，但TP53卻沒(méi)有明顯差異[12].Wang等[13]通過(guò)二代測(cè)序數(shù)據(jù)分析比較了ESCC和EAC樣本在突變、插入缺失、基因擴(kuò)增、純合缺失等基因組改變的異同，發(fā)現(xiàn)在臨床相關(guān)基因組改變中，EAC樣本中KRAS、ERBB2改變更為頻繁，而ESCC樣本中PIK3CA、PTEN、NOTCH1等改變頻率更高.盡管如此，ESCC和EAC中也存在相同的分子變化，如CDKN2A、EGFR、KRAS、MYC、CDK6和MET等基因位點(diǎn)擴(kuò)增同等的發(fā)生于ESCC和EAC.

基于表達(dá)模式相近的基因具有相似的生物學(xué)功能或相同的調(diào)節(jié)機(jī)制的觀點(diǎn)[14]，WGCNA利用分子間加權(quán)表達(dá)相關(guān)性，衡量不同分子的表達(dá)模式關(guān)系.由此，一方面可以近似還原分子相互作用的無(wú)標(biāo)度網(wǎng)絡(luò)分布關(guān)系，另一方面將復(fù)雜的分子表達(dá)數(shù)據(jù)簡(jiǎn)化為若干功能模塊，以挖掘表達(dá)模式相似及具有相同生物學(xué)功能的分子.WGCNA通過(guò)構(gòu)建分子共表達(dá)模塊，常用于鑒定與疾病進(jìn)展、復(fù)發(fā)、預(yù)后相關(guān)的分子功能模塊[15].應(yīng)用WGCNA分析TCGA數(shù)據(jù)庫(kù)中平滑肌瘤測(cè)序數(shù)據(jù)[16]，構(gòu)建24個(gè)功能模塊，發(fā)現(xiàn)其中的Dark red模塊與疾病復(fù)發(fā)呈正相關(guān)(cor=0.32,P=9e-4)，進(jìn)一步用Cytoscape軟件插件CytoHubba鑒定該模塊的12個(gè)節(jié)點(diǎn)基因，結(jié)合平滑肌瘤患者的臨床生存數(shù)據(jù)，最終確定12個(gè)基因中的3個(gè)基因(CDK4、CCT2和MGAT1)是平滑肌瘤的特異性表達(dá)基因，并應(yīng)用ONCOMINE數(shù)據(jù)庫(kù)中的平滑肌瘤數(shù)據(jù)進(jìn)行了外部驗(yàn)證.同樣，利用WGCNA對(duì)16 383個(gè)肝細(xì)胞癌差異表達(dá)基因構(gòu)建25個(gè)模塊[17]，并分析這25個(gè)模塊與肝細(xì)胞癌(HCC)的癌變過(guò)程(正常肝、無(wú)HCC的肝硬化、肝硬化、低度發(fā)育異常、高級(jí)別發(fā)育異常、極早期和早期、晚期HCC和極晚期HCC)的相關(guān)性，發(fā)現(xiàn)5個(gè)節(jié)點(diǎn)基因(GINS1、TOP2A、BUB1B、ACADM和ARPC4)是肝細(xì)胞癌的預(yù)后分子標(biāo)志物.

本研究基于EC的mRNA表達(dá)譜數(shù)據(jù)，應(yīng)用WGCNA構(gòu)建了19個(gè)模塊，其中的Blue和Turquoise模塊分別與EAC、ESCC具有顯著的相關(guān)性.Blue模塊中的基因在EAC中表達(dá)較高，而Turquoise模塊中的基因在ESCC中明顯高表達(dá).此外還發(fā)現(xiàn)Blue模塊與BMI呈正相關(guān)(cor=0.52,P<0.001)，這與之前研究結(jié)果相一致[18].由于BMI較高的肥胖人群導(dǎo)致胃食管反流，使罹患EAC風(fēng)險(xiǎn)增加，表明WGCNA鑒定的Blue模塊及組成基因是EAC和ESCC具有不同發(fā)病機(jī)制的分子基礎(chǔ).

通過(guò)MM與GS之和確定節(jié)點(diǎn)基因的重要性，本研究發(fā)現(xiàn)EPS8L3和SOX15分別是EAC和ESCC的代表性基因，并應(yīng)用3個(gè)不同的數(shù)據(jù)集驗(yàn)證其準(zhǔn)確性.不僅如此，EPS8L3和SOX15還能區(qū)分LUSC和LUAD，進(jìn)一步表明這兩個(gè)分子具有不同的組織起源，是鑒別鱗狀上皮或腺上皮來(lái)源腫瘤的特異性分子，為鱗癌和腺癌的臨床治療方案選擇提供理論基礎(chǔ).EPS8L3與多種腫瘤發(fā)生發(fā)展、化療敏感性相關(guān)，如肝癌組織和細(xì)胞系中EPS8L3高表達(dá)并與肝癌患者的預(yù)后不良有顯著相關(guān)性[19](P<0.011)，EPL8L3基因沉默能夠增加對(duì)順鉑的敏感性.SOX15是一種必要轉(zhuǎn)錄因子，高表達(dá)于LUSC、睪丸生殖細(xì)胞腫瘤、胸腺瘤等，主要通過(guò)調(diào)控WNT/β-catenin信號(hào)通路參與多種腫瘤的發(fā)生、發(fā)展[20].SOX15基因沉默能夠降低WNT/β-catenin信號(hào)通路活性，抑制ESCC細(xì)胞增殖、增加ESCC對(duì)紫杉醇的敏感性[21].上述結(jié)果表明，EPS8L3和SOX15不僅能夠甄別ESCC和EAC，還是潛在分子治療靶點(diǎn).

本研究應(yīng)用WGCNA構(gòu)建ESCC和EAC差異表達(dá)基因的共表達(dá)網(wǎng)絡(luò)，鑒定ESCC和EAC特異表達(dá)基因SOX15和EPS8L3，能夠鑒別起源于食管或肺組織的腫瘤組織類型，并有助于ESCC和EAC的臨床治療方案制定.