通過權(quán)重共表達(dá)網(wǎng)絡(luò)分析胰腺癌轉(zhuǎn)移的特異基因網(wǎng)絡(luò)

唐 婧，程美玲，肖玉鳳，鄧志芬，李 菁

(貴州師范大學(xué)，貴州省植物生理與發(fā)育重點(diǎn)實驗室，生命科學(xué)學(xué)院，貴州 貴陽 550025)

0 引言

胰腺癌(Pancreatic cancer)是一種高度惡性的消化道腫瘤，早期診斷極為困難，預(yù)后生存率極差(5年僅約5%)，病死率近100%，是“癌中之王”[1]。近年來，胰腺癌發(fā)病率逐年升高，癌癥致死已位列第5[2]，預(yù)計到2030年，其將成為第二大致死腫瘤。臨床統(tǒng)計，85%的胰腺癌患者發(fā)生轉(zhuǎn)移，失去手術(shù)治療機(jī)會，甚至有研究認(rèn)為腫瘤形成之前胰腺癌可能已經(jīng)發(fā)生轉(zhuǎn)移[3]，胰腺癌轉(zhuǎn)移調(diào)控機(jī)制、相關(guān)標(biāo)志物及靶向治療是目前胰腺癌研究的重點(diǎn)熱點(diǎn)之一。其中，淋巴結(jié)轉(zhuǎn)移是胰腺癌的重要轉(zhuǎn)移方式，且是影響胰腺癌患者預(yù)后的極其重要因素[4-5]。同時發(fā)現(xiàn)，淋巴結(jié)轉(zhuǎn)移是影響胰腺癌預(yù)后的獨(dú)立因素[6]，無淋巴結(jié)轉(zhuǎn)移者生存率顯著高于有淋巴結(jié)轉(zhuǎn)移患者[7-8]。因此，探索分析胰腺癌轉(zhuǎn)移的特異表達(dá)基因，挖掘相關(guān)標(biāo)志基因及靶標(biāo)基因?qū)τ谝认侔┑念A(yù)后治療具有不可替代且十分重要的意義及作用。

2006年，美國National Cancer Institute(NCI)和National Human Genome Research Institute(NHGRI)聯(lián)合啟動癌癥基因組圖譜項目(The Cancer Genome Atlas，TCGA)，網(wǎng)址：https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga。截止目前為止，該項目數(shù)據(jù)庫收錄了約20 000個癌癥病人，含33類癌癥的約2.5P的基因組、表觀基因組、轉(zhuǎn)錄組和蛋白質(zhì)組數(shù)據(jù)。TCGA數(shù)據(jù)庫中包括177個胰腺癌轉(zhuǎn)錄組數(shù)據(jù)，其包含對應(yīng)的臨床數(shù)據(jù)、測序讀長、拷貝數(shù)變異和生物樣本，及絕大部分的簡單核苷酸變異、DNA甲基化和蛋白組學(xué)數(shù)據(jù)。該數(shù)據(jù)庫為我們分析胰腺癌的調(diào)控機(jī)制、相關(guān)標(biāo)志物及靶向治療提供了豐富可靠的數(shù)據(jù)基礎(chǔ)。

人類基因組計劃后，高通量測序獲取大量癌癥基因組數(shù)據(jù)，通過全基因組關(guān)聯(lián)分析(Genome-wide association study，GWAS)研究使得我們在全基因組范圍內(nèi)獲取大量單個基因與癌癥之間關(guān)系的信息，成功為癌癥的基因診斷，個體化醫(yī)療及靶向治療帶來了契機(jī)。但是GWAS主要基于大樣本從單個基因出發(fā)分析差異基因，而忽略了生物系統(tǒng)的功能作用，缺乏從多基因協(xié)同角度破解癌癥的發(fā)生與發(fā)展。權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析(Weighted gene co-expression network analysis，WGCNA)基于變量相關(guān)性，挖掘生物體系內(nèi)基因互作網(wǎng)絡(luò)體系，利用基因芯片、轉(zhuǎn)錄組數(shù)據(jù)尋找具有協(xié)同表達(dá)特性的基因模塊，進(jìn)一步可以探索基因網(wǎng)絡(luò)與研究表型之間的特異關(guān)系[9]。因此，本研究基于TCGA中胰腺癌數(shù)據(jù)，利用WGCNA挖掘胰腺癌轉(zhuǎn)移相關(guān)的基因網(wǎng)絡(luò)，發(fā)現(xiàn)胰腺癌轉(zhuǎn)移中特異表達(dá)的關(guān)鍵基因，進(jìn)一步分析和胰腺癌轉(zhuǎn)移相關(guān)的代謝通路及蛋白互作，為胰腺癌轉(zhuǎn)移的發(fā)生發(fā)展機(jī)制建立相關(guān)的基因功能和信號通路網(wǎng)絡(luò)分析，以期獲得判斷胰腺癌轉(zhuǎn)移預(yù)后的潛在靶標(biāo)。

1 材料與方法

1.1 數(shù)據(jù)來源

本研究共下載TCGA數(shù)據(jù)庫中的177個胰腺癌組織樣本轉(zhuǎn)錄組數(shù)據(jù)，去除不確定是否轉(zhuǎn)移的樣本后共169個樣本數(shù)據(jù)，其中120個轉(zhuǎn)移樣本(N1)，49個未轉(zhuǎn)移樣本(N0)。

1.2 篩選胰腺癌轉(zhuǎn)移特異差異基因

利用R的DESeq2[10]包進(jìn)行差異表達(dá)基因篩選，用LRT檢驗，篩選padj<0.05的差異基因。

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

利用R的WGCNA[11]算法構(gòu)建胰腺癌權(quán)重共表達(dá)網(wǎng)絡(luò)[9,12]。首先需要構(gòu)建接近真實生物網(wǎng)絡(luò)狀態(tài)的無尺度網(wǎng)絡(luò)，根據(jù)R2>0.85，使用函數(shù)pickSoftThreshold選擇合適的閾值參數(shù)β。然后根據(jù)確定閾值，根據(jù)基因間皮爾遜相關(guān)系數(shù)構(gòu)建鄰接矩陣，轉(zhuǎn)換拓?fù)渚仃嚕又M(jìn)行層次聚類，再用分支切割法(Branch cutting method)[13]來識別構(gòu)建模塊。

1.4 特異基因網(wǎng)絡(luò)模塊的分析

WGCNA算法是以基因網(wǎng)絡(luò)模塊為單位分析特異基因，比較確認(rèn)模塊的保守性，并以此為依據(jù)分析基因網(wǎng)絡(luò)模塊與胰腺癌轉(zhuǎn)移特征之間的關(guān)系，挖掘胰腺癌轉(zhuǎn)移特異的基因網(wǎng)絡(luò)、代謝通路。本研究中利用交叉分析法和保守性統(tǒng)計量計算法比較模塊保守性，分析出胰腺癌轉(zhuǎn)移的特異基因網(wǎng)絡(luò)模塊[14]。交叉分析法是通過Fisher精確檢驗兩個模塊的重疊基因數(shù)，根據(jù)顯著性p值判斷兩個模塊的保守程度。保守性統(tǒng)計量計算法是通過疊加多個與網(wǎng)絡(luò)性質(zhì)相關(guān)的統(tǒng)計量得出一個綜合的統(tǒng)計量Zsummary和MedianRank[14]分析模塊得到保守性。一般Zsummary≥10說明兩個模塊高度保守，2

1.5 模塊基因分析

計算模塊內(nèi)基因的MM和GS值，|MM|>0.9且|GS|>0.2的基因為網(wǎng)絡(luò)模塊的樞紐基因，它們對具有模塊主要的調(diào)控作用。用Cytoscape軟件繪制模塊中共表達(dá)基因網(wǎng)絡(luò)關(guān)系[15]。利用超幾何分布的修正Fisher精確檢驗基因在KEGG上的富集[16]。

2 結(jié)果與分析

2.1 表達(dá)差異基因

胰腺癌轉(zhuǎn)移轉(zhuǎn)錄組樣本120個(N1)，未發(fā)生轉(zhuǎn)移胰腺癌轉(zhuǎn)錄組樣本49個(N0)，共有55 126個基因，通過DESeq篩選得到1 590個表達(dá)差異基因(padj<0.05，log2FC(N1/N0)>2，F(xiàn)C：FoldChange)。

圖1 差異表達(dá)基因火山圖Fig.1 The volcano plot of differentially expressed genes

2.2 WGCNA分析

2.2.1 網(wǎng)絡(luò)構(gòu)建及模塊識別

首先，對胰腺癌轉(zhuǎn)移樣本(N1)和未轉(zhuǎn)移樣本(N0)利用Hclust函數(shù)進(jìn)行樣品聚類分析，根據(jù)閾值分別去除N1和N0中的7個離群樣本；然后根據(jù)R2>0.85的要求，選取的β值(β=12)；根據(jù)β使用Adjacency函數(shù)生成鄰接矩陣，TOMsimilarity函數(shù)將鄰接矩陣轉(zhuǎn)化為拓?fù)渲丿B矩陣，使用動態(tài)剪枝圖法進(jìn)行動態(tài)分支切割，設(shè)置每個模塊中最少包含基因數(shù)為30，構(gòu)建網(wǎng)絡(luò)模塊。最終N1和N0分別構(gòu)建為6個模塊和8個模塊(圖2)，其中，Grey模塊代表未被分到其他任何模塊的基因。胰腺癌轉(zhuǎn)移特異差異基因分別在轉(zhuǎn)移(N1)和未轉(zhuǎn)移(N0)樣本中，構(gòu)建成6和8個基因共表達(dá)網(wǎng)絡(luò)模塊，說明在胰腺癌轉(zhuǎn)移的過程中，極大的加強(qiáng)了基因之間的共表達(dá)特性，增強(qiáng)了基因間的互作作用，基因共活躍程度增高，相互影響作用加大。

圖2 基因聚類樹和模塊劃分(左N1，右N0)Fig.2 Clustering dendrogram of genes and module division (Left: N1, Right: N0)

2.2.2 模塊比較分析

通過交叉比較法，比較分析胰腺癌轉(zhuǎn)移與未轉(zhuǎn)移之間的差異，結(jié)果見3。N0組網(wǎng)絡(luò)中識別出8個模塊，而N1組中只識別出了6個。首先，胰腺癌轉(zhuǎn)移后，合并模塊Grey模塊內(nèi)的基因數(shù)(9)明顯少于未轉(zhuǎn)移胰腺癌(162)，幾乎全部胰腺癌轉(zhuǎn)移差異表達(dá)基因全部高相關(guān)度的構(gòu)建到基因共表達(dá)模塊中；說明在胰腺癌轉(zhuǎn)移過程中，差異表達(dá)基因的相互影響作用顯著性高于未轉(zhuǎn)移胰腺癌。然后，胰腺癌樣本N0中的Black、Magenta和Pink模塊在胰腺癌轉(zhuǎn)移樣本N1中全部并入最大Turquoise基因共表達(dá)網(wǎng)絡(luò)模塊(1 013個基因)；進(jìn)一步表明在胰腺癌轉(zhuǎn)移后基因相互作用加大，基因間共表達(dá)效應(yīng)更強(qiáng)，合并形成更大的共表達(dá)互作網(wǎng)絡(luò)。除此之外，胰腺癌轉(zhuǎn)移樣本中的Yellow網(wǎng)絡(luò)模塊與未轉(zhuǎn)移樣本中的網(wǎng)絡(luò)模塊幾乎不存在顯著性基因重疊，說明在胰腺癌轉(zhuǎn)移過程中，基因在互作加強(qiáng)的基礎(chǔ)上，構(gòu)建出一個新的基因共表達(dá)網(wǎng)絡(luò)關(guān)系；該網(wǎng)絡(luò)內(nèi)含有52個基因。

圖3 轉(zhuǎn)移(N1)模塊與未轉(zhuǎn)移(N0)模塊的交叉法比較(交疊數(shù)值上方數(shù)值表示重疊基因數(shù)目；下方數(shù)值表示p值； 橫坐標(biāo)：N0模塊；縱坐標(biāo)：N1模塊)Fig.3 The cross-tabulation of metastasis modules (columns) and without metastases modules (rows)(The number above the overlap value indicates the number of overlapping genes; the down indicates p value. Abscissa: N0 modules; Ordinate; N1 modules)

進(jìn)一步通過模塊保守性分析，我們比較胰腺癌轉(zhuǎn)移與未轉(zhuǎn)移之間基因網(wǎng)絡(luò)模塊之間的差異，結(jié)果如圖4所示。由圖4可以發(fā)現(xiàn)：去除未分類的Grey模塊，只有Yellow模塊的Z_summary值小于2，說明胰腺癌轉(zhuǎn)移樣本的Yellow基因網(wǎng)絡(luò)模塊與未轉(zhuǎn)移樣本的基因模塊顯著不同，是胰腺癌轉(zhuǎn)移中形成的新的基因共表達(dá)關(guān)系。同時Yellow基因網(wǎng)絡(luò)模塊的MedianRank最高，進(jìn)一步證實了其保守性最低，顯著不同于未轉(zhuǎn)移樣本的基因模塊。

圖4 轉(zhuǎn)移(N1)模塊與未轉(zhuǎn)移(N0)模塊的網(wǎng)絡(luò)的保守性分析Fig.4 The preservation test of metastasis modules (N1) and without metastases modules (N0)

綜上所述，我們推測相對于胰腺癌未轉(zhuǎn)移的樣本，在胰腺癌轉(zhuǎn)移后胰腺癌轉(zhuǎn)移特異表達(dá)基因形成一個新的強(qiáng)共表達(dá)互作關(guān)系的基因網(wǎng)絡(luò)，該網(wǎng)絡(luò)與胰腺癌轉(zhuǎn)移密切有關(guān)；網(wǎng)絡(luò)內(nèi)共包含52個基因。

2.3 Yellow模塊分析

2.3.1 Yellow網(wǎng)絡(luò)及樞紐基因分析

Yellow模塊是由52個胰腺癌轉(zhuǎn)移特意表達(dá)基因構(gòu)成共表達(dá)網(wǎng)絡(luò)(圖5)。樞紐基因是指在模塊中連接度影響值最高的一些基因[17]。結(jié)果發(fā)現(xiàn)，在與胰腺癌轉(zhuǎn)移密切相關(guān)的特有共表達(dá)基因網(wǎng)絡(luò)Iellow中，根據(jù)|MM|和|GS|值篩選，我們發(fā)現(xiàn)Yellow網(wǎng)絡(luò)模塊含有ASIC4，GLDC，MMD2，CTNNA2，F(xiàn)OXC4-AS1，SNHG19和BCAN等11個樞紐基因；它們與胰腺癌轉(zhuǎn)移密切相關(guān)，極可能是潛在的胰腺癌轉(zhuǎn)移治療和診斷的重要靶標(biāo)。

圖5 Yellow模塊基因網(wǎng)絡(luò)及其樞紐基因(紅色外圈代表樞紐基因)Fig.5 The co-expression network of yellow module and hub genes (Red outer circle represents the hub gene)

2.3.2 KEGG通路富集分析

對Yellow模塊中52個基因進(jìn)行KEGG通路富集，結(jié)果見圖6。我們發(fā)現(xiàn)顯著富集于炎癥介質(zhì)對色氨酸通道的調(diào)節(jié)(Inflammatory mediator regulation of TRP channels)，卵母細(xì)胞減數(shù)分裂(Oocyte meiosis)，乙醛酸和二元羧酸鹽代謝(Glyoxylate and dicarboxylate metabolism)和甘氨酸、絲氨酸和蘇氨酸代謝(Glycine, serine and threonine metabolism)等4條代謝通路，推測認(rèn)為它們可能與胰腺癌轉(zhuǎn)移密切相關(guān)。

圖6 KEGG通路富集分布情況Fig.6 KEGG pathway enrichment

3 討論

近幾年來，隨著高通量測序技術(shù)的普及，越來越多的復(fù)雜疾病基礎(chǔ)研究開始應(yīng)用轉(zhuǎn)錄組數(shù)據(jù)來分析潛在的致病基因。除了應(yīng)用最為廣泛的差異表達(dá)基因分析之外，依據(jù)基因與基因之間的關(guān)聯(lián)關(guān)系的共表達(dá)網(wǎng)絡(luò)分析也逐漸被許多研究者們采納。本項目對TCGA上的胰腺癌轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行權(quán)重基因共表達(dá)網(wǎng)絡(luò)分析，成功定位出一個與胰腺癌轉(zhuǎn)移密切相關(guān)的新基因共表達(dá)網(wǎng)絡(luò)模塊，模塊內(nèi)含52個胰腺癌差異表達(dá)基因；且找到ASIC4，GLDC，MMD2，CTNNA2，F(xiàn)OXC4-AS1，SNHG19和BCAN等11個樞紐基因調(diào)控該基因網(wǎng)絡(luò)模塊，它們極可能與胰腺癌轉(zhuǎn)移調(diào)控有重要關(guān)系及作用。

在過往的研究中已驗證CTNNA2為新型miR-4885的靶基因，它在食管癌組織中的表達(dá)上調(diào)[17]。CTNNA2兩個SNP與印尼女性患乳腺癌的風(fēng)險呈暗示性關(guān)聯(lián)[18]，且與胰腺癌存活之間可能存在機(jī)械聯(lián)系[19]。FOXC2-AS1已被確認(rèn)為骨肉瘤中的癌基因[20]。還發(fā)現(xiàn)模塊中SMC1B與胰腺癌中最常見的腫瘤抑制因子p53的通路網(wǎng)絡(luò)連接，并與胰腺疾病建立聯(lián)系[21]。RGS 5、BCAN、MAP2K6可能被用作這些癌癥中的新型診斷或預(yù)后標(biāo)志物[22-24]。

綜上所述，本項目通過WCGNA方法利用胰腺癌轉(zhuǎn)移和非轉(zhuǎn)移樣本之間的基因關(guān)聯(lián)關(guān)系的差異預(yù)測，有效的分析了與胰腺癌轉(zhuǎn)移有關(guān)的潛在基因模塊。