基于加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析構(gòu)建胃癌轉(zhuǎn)移預(yù)測(cè)模型

龔 超，陳 魁，章德昆，謝徑峰，吳芳華，黃玉鈿，薛玉欽，王力群

1.福建醫(yī)科大學(xué)附屬福州市第一醫(yī)院普通外科，福建 福州 350009；2.福建醫(yī)科大學(xué)附屬福州市第一醫(yī)院檢驗(yàn)科，福建 福州 350009；3.福建醫(yī)科大學(xué)附屬福州市第一醫(yī)院病理科，福建 福州 350009

胃癌是世界上癌癥相關(guān)性死亡的主要原因之一，中國(guó)2015年胃癌發(fā)病率和死亡率中均居癌癥發(fā)病率和相關(guān)死亡率的第二位［1］。盡管目前手術(shù)結(jié)合化療的方案已廣泛應(yīng)用于胃癌的治療當(dāng)中，但胃癌的5年生存率依然較低［2-3］。據(jù)文獻(xiàn)［4-5］報(bào)道，非轉(zhuǎn)移性胃癌患者5年生存率可達(dá)60%，但晚期轉(zhuǎn)移性胃癌患者5年生存率不足10%。因此，早期識(shí)別具有高危轉(zhuǎn)移因素的胃癌患者，采取更積極的治療方案具有重要的臨床價(jià)值。

近年來(lái)基因芯片篩選技術(shù)廣泛應(yīng)用于臨床治療，同時(shí)有研究［6-7］表明，利用基因芯片技術(shù)篩選胃癌轉(zhuǎn)移相關(guān)基因具有重要價(jià)值。同時(shí)每例患者的生物調(diào)控機(jī)制涉及復(fù)雜的基因網(wǎng)絡(luò)互作，僅通過(guò)傳統(tǒng)的基因芯片篩選差異基因不可否認(rèn)其可對(duì)生物系統(tǒng)的局部作出解釋，但難免遺漏調(diào)控過(guò)程中的核心分子。本研究利用高通量功能基因組數(shù)據(jù)庫(kù)（Gene Expression Omnibus，GEO）中含有晚期轉(zhuǎn)移胃癌數(shù)據(jù)的數(shù)據(jù)集（GES103236），采用加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析（weighted gene co-expression network analysis，WGCNA）方法篩選出與胃癌轉(zhuǎn)移相關(guān)的靶基因并構(gòu)建預(yù)測(cè)模型，擬為胃癌患者個(gè)體化診療提供參考依據(jù)。

1 資料和方法

1.1 GEO數(shù)據(jù)集獲取

本研究從GEO（https://www.ncbi.nlm.nih.gov/geo/）中，獲取了微陣列數(shù)據(jù)集（GSE103236和GSE14210），含有10例轉(zhuǎn)移和非轉(zhuǎn)移的胃癌組織，以及9例癌旁組織，我們使用該數(shù)據(jù)集作為建模數(shù)據(jù)集，用于構(gòu)建基因共表達(dá)網(wǎng)絡(luò)模型。GSE14210含有的123例有疾病進(jìn)展和預(yù)后情況患者用于建立預(yù)測(cè)模型驗(yàn)證患者疾病進(jìn)展情況。我們利用Kaplan-Meier plot網(wǎng)站（http://kmplot.com/analysis/）分析多個(gè)GEO數(shù)據(jù)集中胃癌患者的生存數(shù)據(jù)。同時(shí)，含有多種癌癥組織和癌旁組織的ONCOMINE數(shù)據(jù)庫(kù)（https://www.oncomine.org）也被用于進(jìn)行外部驗(yàn)證。

1.2 WGCNA

WGCNA是一種常用的模塊化分析技術(shù)，已被用于識(shí)別和篩選復(fù)雜疾病的生物標(biāo)志物或藥物靶點(diǎn)［8］。首先，我們對(duì)樣本的基因表達(dá)情況進(jìn)行質(zhì)檢以及離群值分析，以確保均為樣本且能夠正常使用。進(jìn)一步，我們通過(guò)R軟件中的“WGCNA”分析包構(gòu)建基因共表達(dá)網(wǎng)絡(luò)［9-10］。然后，通過(guò)主要連接關(guān)系和皮爾森相關(guān)矩陣建立兩個(gè)基因之間的相關(guān)矩陣。并通過(guò)對(duì)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的分析，確定軟閾值大小。軟閾值是基于近似無(wú)標(biāo)度拓?fù)涞囊环N準(zhǔn)則。無(wú)標(biāo)度拓?fù)鋽M合指數(shù)曲線在達(dá)到較高值后趨于平緩，選擇的軟閾值還需要有一個(gè)更好的平均連通性。我們進(jìn)一步將鄰接轉(zhuǎn)化為拓?fù)渲丿B矩陣（topological overlap matrix，TOM），TOM可以度量一個(gè)基因的網(wǎng)絡(luò)連通性，該基因定義為其與其他所有基因鄰接的總和，用于網(wǎng)絡(luò)生成［11-12］。為了將表達(dá)譜相似的基因分類到基因模塊中，根據(jù)基于TOM的差異測(cè)度進(jìn)行平均連鎖層次聚類［9-10］。模塊鑒定后，根據(jù)各組表型數(shù)據(jù)，采用t檢驗(yàn)計(jì)算各組間各基因表達(dá)顯著性檢驗(yàn)的P值。每個(gè)模塊的顯著性定義為模塊內(nèi)基因顯著性的平均值。具有顯著性的模塊可能與特定疾病的存在有關(guān)。為了進(jìn)一步分析模塊，我們計(jì)算模塊特征基因（MEs）的差異性，為模塊樹狀圖選擇一條切線，并合并部分模塊。

1.3 特征模塊和核心基因篩選

在對(duì)每個(gè)基因模塊進(jìn)行主分析時(shí)，MEs被視為主要成分，所有基因的表達(dá)模式可歸納為給定模塊內(nèi)的單一特征表達(dá)譜。另外，我們通過(guò)皮爾森相關(guān)檢驗(yàn)來(lái)評(píng)估MEs與轉(zhuǎn)移的相關(guān)性，以確定相關(guān)模塊。選擇與轉(zhuǎn)移高度相關(guān)的模塊作為轉(zhuǎn)移模塊進(jìn)一步進(jìn)行分析。此外，為了確定與轉(zhuǎn)移相關(guān)的模塊，我們進(jìn)行了模塊與轉(zhuǎn)移之間的皮爾森相關(guān)性分析。識(shí)別中心的基因，我們選擇轉(zhuǎn)移模型相關(guān)系數(shù)最高的數(shù)據(jù)集，該模塊也是在所有的模塊中比重最大，在模塊和中心基因定義的模塊連接以絕對(duì)值衡量皮爾森的相關(guān)性。為了尋找真正的核心靶基因，我們對(duì)數(shù)據(jù)集進(jìn)行了差異基因分析，同時(shí)利用VEEN圖對(duì)癌組織和癌旁組織的差異基因、轉(zhuǎn)移和非轉(zhuǎn)移的差異基因和Lightsteelblue模塊中基因取交集。我們篩選出了C5AR1、AP3M2、TYMP、ANXA2P1用于進(jìn)一步分析。

1.4 基因本體（gene ontology，GO）功能學(xué)和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）通路分析

為了解Lightsteelblue模塊主要涉及的功能學(xué)和通路。我們采用標(biāo)準(zhǔn)富集計(jì)算方法進(jìn)行GO功能分析和KEGG通路分析用以篩選與其相關(guān)的功能和通路。

1.5 統(tǒng)計(jì)學(xué)處理

統(tǒng)計(jì)分析采用的是SPSS 24.0，GraphPad Prism 7.0和R 3.4.1，同時(shí)用上述軟件進(jìn)行圖像生成處理。t檢驗(yàn)用來(lái)分析兩個(gè)組別之間的平均數(shù)的差異。繪制受試者工作特征（receiver operating characteristic，ROC）曲線，評(píng)估核心靶基因的預(yù)測(cè)能力，利用曲線下面積（area under curve，AUC）評(píng)估靈敏度和特異度。利用Kaplan-Meier法繪制生存曲線，預(yù)測(cè)基因?qū)Σ』碱A(yù)后的影響。P<0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 差異基因分析

利用“edge”R篩選差異表達(dá)的基因，基于fold change=1（P<0.05）為閾值篩選GSE103236數(shù)據(jù)集中癌組織和癌旁組織差異基因，共發(fā)現(xiàn)3 431個(gè)差異表達(dá)基因，其中2 399個(gè)基因表達(dá)上調(diào)，1 032個(gè)基因表達(dá)下調(diào)。同理，對(duì)轉(zhuǎn)移和非轉(zhuǎn)移患者進(jìn)行差異基因篩選，共發(fā)現(xiàn)1 264個(gè)差異表達(dá)基因，其中1 218個(gè)基因表達(dá)上調(diào)，46個(gè)基因表達(dá)下調(diào)。同時(shí)，通過(guò)R軟件繪制熱圖，50個(gè)基因在10個(gè)癌組織、9個(gè)癌旁組織、7個(gè)非轉(zhuǎn)移癌組織和3個(gè)轉(zhuǎn)移癌組織中的表達(dá)情況見(jiàn)圖1。

圖1 熱圖分析 Fig.1 Thermographic analysis

2.2 構(gòu)建共表達(dá)模塊

基于19 711個(gè)差異基因在10例胃癌組織和9例癌旁組織的表達(dá)數(shù)據(jù)，利用基于無(wú)序列網(wǎng)絡(luò)的WGCNA方法將基因進(jìn)行模塊化富集分析（圖2），將基因依據(jù)其各相關(guān)表達(dá)量進(jìn)一步進(jìn)行分類。為此，共篩選獲得了79個(gè)相應(yīng)的基因模塊，為了進(jìn)一步了解各模塊與轉(zhuǎn)移的關(guān)聯(lián)度，本研究進(jìn)一步將臨床信息納入分析，檢測(cè)各模塊與轉(zhuǎn)移的皮爾森相關(guān)系數(shù)，以便篩選出與轉(zhuǎn)移最相關(guān)的模塊，納入進(jìn)一步的分析。我們依據(jù)各模塊在轉(zhuǎn)移上皮爾森系數(shù)絕對(duì)值相加為最高者認(rèn)定為響應(yīng)系數(shù)最高模塊，最后篩選獲取了Lightsteelblue模塊。

圖2 WGCNA模式圖Fig.2 WGCNA cluster dendrogram

2.3 基因模塊與臨床性狀關(guān)聯(lián)分析

為尋找與轉(zhuǎn)移相關(guān)的關(guān)鍵調(diào)控基因及效應(yīng)通路，基于GO功能富集以及KEGG數(shù)據(jù)庫(kù)中信號(hào)轉(zhuǎn)導(dǎo)通路的上下游關(guān)系，本課題組依據(jù)Lightsteelblue模塊中所含有的47個(gè)相關(guān)基因篩選表達(dá)響應(yīng)基因參與的信號(hào)通路的關(guān)聯(lián)圖（圖3A、B）。我們據(jù)此發(fā)現(xiàn)多條極為相關(guān)的信號(hào)轉(zhuǎn)導(dǎo)通路，如細(xì)胞黏附、轉(zhuǎn)錄異常調(diào)控等。目前已有較多的研究［13-15］顯示，上述通路在胃癌中涉及轉(zhuǎn)移敏感性。因此初步證明了本課題前期中所篩選出的Lightsteelblue模塊是和轉(zhuǎn)移相關(guān)的模塊。同時(shí)，為了進(jìn)一步驗(yàn)證該模塊是否與轉(zhuǎn)移相關(guān)，我們將Lightsteelblue模塊中基因與轉(zhuǎn)移一起予以分析（圖4），結(jié)果顯示，在Lightsteelblue模塊中大部分基因其P值均小于0.05，進(jìn)一步佐證了Lightsteelblue模塊與轉(zhuǎn)移的相關(guān)性?；诖宋覀兛梢猿醪秸J(rèn)定Lightsteelblue模塊就是與轉(zhuǎn)移關(guān)聯(lián)最大的模塊。為進(jìn)一步確保模塊中的基因與胃癌發(fā)生、發(fā)展及轉(zhuǎn)移相關(guān)，我們進(jìn)一步通過(guò)VEEN分析，對(duì)模塊中的基因和我們前期所篩選的差異基因取交集。最終，我們共篩選了4個(gè)基因C5AR1、AP3M2、TYMP、ANXA2P1作為核心靶基因（圖5）。

圖3 GO分析和KEGG分析Fig.3 GO analysis and KEGG analysis

圖4 Lightsteelblue模塊基因與轉(zhuǎn)移相關(guān)分析Fig.4 Analysis of the relationship between Lightsteelblue module gene and metastasis

圖5 VEEN圖Fig.5 VEEN

2.4 核心靶基因驗(yàn)證

為了進(jìn)一步驗(yàn)證我們篩選的基因?qū)︻A(yù)測(cè)患者胃癌是否轉(zhuǎn)移以及預(yù)后的效能，我們重新納入原始數(shù)據(jù)集進(jìn)一步分析。本課題前期篩選的4個(gè)基因表達(dá)在腫瘤組織中表達(dá)顯著高于癌旁組織（P<0.01，圖6），同時(shí)在轉(zhuǎn)移組中，其表達(dá)量顯著高于非轉(zhuǎn)移組。進(jìn)一步證明了我們前期篩選的基因不僅與胃癌轉(zhuǎn)移相關(guān)，而且為明確的癌基因。為驗(yàn)證基因表達(dá)豐度能否判定腫瘤的發(fā)生或者轉(zhuǎn)移，采用了ROC曲線計(jì)算其AUC用于判定各基因在預(yù)測(cè)腫瘤發(fā)生以及預(yù)測(cè)轉(zhuǎn)移的概率。

圖6 各基因相對(duì)表達(dá)量Fig.6 The relative expressions of different genes

在預(yù)測(cè)腫瘤的發(fā)生中，C5AR1、AP3M2、TYMP、ANXA2P1基因其分別預(yù)測(cè)概率為0.767、0.844、0.956和0.867?；诖宋覀冞M(jìn)一步建立了回歸預(yù)測(cè)模型，該模型預(yù)測(cè)腫瘤發(fā)生其AUC可達(dá)1。上述結(jié)果提示前期所篩選的基因在預(yù)測(cè)腫瘤發(fā)生中也有較大的意義。為進(jìn)一步驗(yàn)證其用于預(yù)測(cè)胃癌轉(zhuǎn)移的效能，在預(yù)測(cè)胃癌轉(zhuǎn)移發(fā)生中，C5AR1、AP3M2、TYMP、ANXA2P1基因其分別預(yù)測(cè)概率為0.952、0.905、0.952和0.857。基于此我們進(jìn)一步建立了LOGSTICAL預(yù)測(cè)模型，該模型預(yù)測(cè)腫瘤發(fā)生其AUC可達(dá)1。上述結(jié)果說(shuō)明我們所篩選的基因與胃癌發(fā)生、發(fā)展和轉(zhuǎn)移相關(guān)（圖7）。

圖7 ROC曲線預(yù)測(cè)腫瘤發(fā)生和轉(zhuǎn)移 Fig.7 The prediction of tumor genesis and metastasis by ROC curves

2.5 核心靶基因外部表達(dá)和生存驗(yàn)證

我們利用ONCOMINE數(shù)據(jù)庫(kù)中包含有大量的mRNA測(cè)序/芯片結(jié)果，對(duì)C5AR1、AP3M2、TYMP和ANXA2P1基因多數(shù)據(jù)集的meta分析，在所納入的8個(gè)數(shù)據(jù)集中，所篩選的C5AR1、AP3M2、TYMP和ANXA2P1基因在癌組織中均相對(duì)癌旁組織呈現(xiàn)高表達(dá)（P=0.002，P=0.001，P<0.001，P=0.003，圖8）。結(jié)果表明我們的前期結(jié)果是可信的。同時(shí)我們利用Kaplan-Meier plot網(wǎng)站（http://kmplot.com/analysis/）所含有的多個(gè)胃癌生存數(shù)據(jù)集分析我們所篩選的基因是否可以預(yù)測(cè)預(yù)后。結(jié)果顯示，在胃癌患者中，C5AR1、AP3M2、TYMP和ANXA2P1基因低表達(dá)患者有著更好的預(yù)后，進(jìn)一步佐證了我們前期結(jié)果（圖9）。

圖8 ONCOMINE數(shù)據(jù)庫(kù)分析圖Fig.8 ONCOMINE database analysis

圖9 Kaplan-Meier plot數(shù)據(jù)庫(kù)分析生存情況Fig.9 Survival curves by Kaplan-Meier plot database analysis

2.6 利用GSE14210構(gòu)建核心基因預(yù)測(cè)復(fù)發(fā)和生存模型

為進(jìn)一步明確所篩選的基因是否能有效地預(yù)測(cè)胃癌患者的復(fù)發(fā)及生存情況，我們利用生存風(fēng)險(xiǎn)模型分析GEO數(shù)據(jù)庫(kù)中帶有總生存以及疾病進(jìn)展資料的數(shù)據(jù)庫(kù)GSE14210。其內(nèi)共有167例胃癌患者基因芯片數(shù)據(jù)，有123例患者帶有確切的疾病進(jìn)展以及生存資料數(shù)據(jù)。本研究利用該數(shù)據(jù)集，通過(guò)對(duì)上述核心基因表達(dá)數(shù)據(jù)進(jìn)行提取以及生存數(shù)據(jù)分析，構(gòu)建生存風(fēng)險(xiǎn)預(yù)測(cè)模型（圖10），其公式=C5AR1×（0.113 8）+AP3M2×（0.000 1）+TYMP×（0.093 9）+ANXA2P1×（0.486 3）。據(jù)此，我們可以獲得每個(gè)患者的風(fēng)險(xiǎn)得分，并利用該得分對(duì)患者是否出現(xiàn)疾病進(jìn)展進(jìn)行預(yù)測(cè)，結(jié)果表明，所構(gòu)建的生存風(fēng)險(xiǎn)模型能較好地預(yù)測(cè)胃癌患者治療后是否出現(xiàn)疾病進(jìn)展（AUC=0.71，P<0.001，圖11）。據(jù)此有理由相信我們所篩選的基因以及模型是可靠的。

圖10 構(gòu)建生存風(fēng)險(xiǎn)模型Fig.10 Survival risk model

圖11 ROC曲線Fig.11 ROC curve

3 討 論

胃癌是世界上癌癥相關(guān)性死亡的主要原因之一，盡管目前手術(shù)結(jié)合化療/免疫治療/靶向治療等多種方案已廣泛應(yīng)用于胃癌的治療當(dāng)中，但胃癌的5年生存率依然低下［2-3］。這主要是由于大多數(shù)胃癌被發(fā)現(xiàn)時(shí)已處于中晚期，所以預(yù)后不佳。晚期轉(zhuǎn)移性胃癌患者5年生存率不足 10%［4-5］，因此，本研究為更早識(shí)別具有高危轉(zhuǎn)移風(fēng)險(xiǎn)的胃癌患者，對(duì)GEO數(shù)據(jù)集轉(zhuǎn)移和非轉(zhuǎn)移胃癌組織及癌旁組織的數(shù)據(jù)集，通過(guò)精準(zhǔn)的WGCNA算法識(shí)別出了4種與胃癌發(fā)生、轉(zhuǎn)移相關(guān)的基因，即C5AR1、AP3M2、TYMP和ANXA2P1；同時(shí)通過(guò)內(nèi)部表達(dá)驗(yàn)證和利用外部的ONCOMINE數(shù)據(jù)庫(kù)和Kaplan-Meier plot數(shù)據(jù)庫(kù)驗(yàn)證了我們所篩選的基因具有較強(qiáng)的重復(fù)性和可靠性。

WGCNA通過(guò)將具有相似的功能或者出現(xiàn)在同一生物通路中的基因視作一個(gè)模塊，進(jìn)而在基因模塊水平上研究每個(gè)生物個(gè)體與外部信息的聯(lián)系。并通過(guò)有效地劃分基因模塊進(jìn)一步研究生物調(diào)控網(wǎng)絡(luò)內(nèi)部關(guān)聯(lián)性［8］。因此基于WGCNA的方法構(gòu)建基因網(wǎng)絡(luò)互作關(guān)系模式，不僅可以幫助我們從基因網(wǎng)絡(luò)調(diào)控角度更為精準(zhǔn)地篩選出核心靶基因，又可避免由于患者個(gè)體基因表達(dá)差異導(dǎo)致的假陽(yáng)/陰性結(jié)果。本研究利用該方法，篩選出了C5AR1、AP3M2、TYMP和ANXA2P1基因。同時(shí)，我們通過(guò)表達(dá)分析發(fā)現(xiàn)，這些基因在胃癌組織中相較于癌旁組織呈高表達(dá)，在轉(zhuǎn)移患者中表達(dá)量較非轉(zhuǎn)移患者高。上述結(jié)果說(shuō)明我們通過(guò)WGCNA方法篩選的基因是有效的。為了在外部數(shù)據(jù)庫(kù)中驗(yàn)證，我們利用了ONCOMINE數(shù)據(jù)庫(kù)，其包含多個(gè)胃癌和癌旁組織測(cè)序/基因芯片的結(jié)果，通過(guò)對(duì)多個(gè)數(shù)據(jù)集進(jìn)行的meta分析發(fā)現(xiàn)，上述4個(gè)基因在胃癌組織中相較于癌旁組織均高表達(dá)。同時(shí)，我們還通過(guò)包含多個(gè)不同來(lái)源數(shù)據(jù)集胃癌患者生存數(shù)據(jù)的Kaplan-Meier plot網(wǎng)站，驗(yàn)證了C5AR1、AP3M2、TYMP和ANXA2P1基因與患者預(yù)后的關(guān)系。結(jié)果顯示，上述4個(gè)基因高表達(dá)均預(yù)示著患者有不良預(yù)后。同時(shí)，本研究還利用了GSE14210數(shù)據(jù)集，基于上述的核心靶基因構(gòu)建了疾病進(jìn)展和生存模型，該模型的結(jié)果也能較好地預(yù)測(cè)患者疾病進(jìn)展，從而進(jìn)一步佐證了我們所篩選的核心基因可以有效地預(yù)測(cè)患者的疾病進(jìn)展和預(yù)后。因此，本研究結(jié)果對(duì)未來(lái)闡明C5AR1、AP3M2、TYMP和ANXA2P1基因表達(dá)與胃癌轉(zhuǎn)移及預(yù)后的關(guān)系奠定了一定的基礎(chǔ)。但不可否認(rèn)的是，本研究的結(jié)果還需要在大量臨床實(shí)際樣本中進(jìn)一步驗(yàn)證，并通過(guò)一系列的體內(nèi)/體外實(shí)驗(yàn)探索上述4個(gè)基因的作用、相互關(guān)系、機(jī)制及臨床意義。

總之，C5AR1、AP3M2、TYMP和ANXA2P1基因有可能成為新的預(yù)后指標(biāo)，有助于胃癌患者個(gè)性化治療及臨床預(yù)后判斷。