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

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

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

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

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

1 材料和方法

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

從癌癥基因圖譜數(shù)據(jù)庫(the cancer genome atlas, TCGA)和基因表達數(shù)據(jù)庫(gene expression omnibus, GEO)下載EC的mRNA轉錄組數(shù)據(jù)及相應臨床信息，分別為TCGA-EC和GSE26866.將表達譜數(shù)據(jù)與臨床信息匹配后，共納入來源于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差異表達基因鑒定

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

1.3 WGCNA網(wǎng)絡模塊構建

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

根據(jù)模塊第一主成分Module Eigengenes(ME)計算模塊內所有基因與ME的相關性，確定基因的模塊隸屬度(module membership, MM)；根據(jù)ME與EAC、ESCC表型的相關性，確定基因顯著性(gene significance, GS)；MM與GS均滿足一定條件的基因為模塊的節(jié)點基因．

1.4 富集分析及獨立驗證

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

應用3組獨立數(shù)據(jù)驗證上述關鍵基因鑒別ESCC和EAC的準確性，分別是GSE26866、肺鱗狀細胞癌(lung squamous cell carcinoma, LUSC)和肺腺癌(lung adenocarcinoma, LUAD)、癌細胞數(shù)據(jù)庫(cancer cell line encyclopedia, CCLE)中23個ESCC細胞系和1個EAC細胞系.

1.5 統(tǒng)計學分析

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

2 結果

2.1 EAC和ESCC差異表達基因

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

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

2.2 共表達網(wǎng)絡模塊分析

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

上方為基因聚類圖譜；下方為動態(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 食管鱗癌和食管腺癌差異表達基因和共表達網(wǎng)絡模塊分布Fig.2 Distribution of differentially expressed genes between esophageal squamous cell carcinoma and esophageal adenocarcinoma, and co-expression network modules by WGCNA

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

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

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

2.3 EAC和ESCC代表性基因確定及驗證

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

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

(a-b): blue和turquoise模塊內基因的MM與GS特征散點圖； (c-d): EPS8L3和SOX15在數(shù)據(jù)集TCGA中的表達值比較； (e-f): EPS8L3和SOX15在數(shù)據(jù)集GSE26886中的表達譜比較，差異顯著性為FDR法矯正后的P值； (g-h): EPS8L3和SOX15在數(shù)據(jù)集LUSC、LUAD中的表達譜比較. 兩組樣本比較用Wilcoxon檢驗；多組樣本比較用Kruskal-Wallis檢驗；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 關鍵節(jié)點基因表達值比較Fig.4 Comparation of key hub genes from modules and their expression values.

2.4 基因生物學功能和通路富集分析

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

(a): Blue模塊GO富集分析結果； (b): Turquoise模塊的GO富集分析結果.(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 討論

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

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

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

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

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

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