基于加權(quán)基因共表達(dá)網(wǎng)絡(luò)的膀胱尿路上皮癌預(yù)后及療效相關(guān)基因篩選

李佳麗，曾梓航，劉嘉彬，張 冉，李 碩，孟詳喻，李 勝，李 偉

膀胱尿路上皮癌(bladder urothelial carcinoma，BLCA)占膀胱癌發(fā)病率的90%以上，是泌尿系統(tǒng)中最常見的惡性腫瘤之一[1-2]。BLCA具有多發(fā)性、復(fù)發(fā)率高等特點(diǎn)，其發(fā)病機(jī)制尚未完全明確。目前有研究表明，BLCA的發(fā)病機(jī)制可能和多基因、多層次的細(xì)胞信號網(wǎng)絡(luò)紊亂有關(guān)[3]。BLCA的治療效果(即進(jìn)展、穩(wěn)定、部分緩解、完全緩解)顯著地影響了患者的預(yù)后，且其基因?qū)用娴臋C(jī)制有待闡明。

加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析(weighed gene co-expression network analysis,WGCNA)是一種構(gòu)建基因共表達(dá)網(wǎng)絡(luò)的系統(tǒng)生物學(xué)算法，該算法以轉(zhuǎn)錄組數(shù)據(jù)為分析基礎(chǔ)，在世界生物醫(yī)學(xué)范圍均得到廣泛的應(yīng)用[4]。WGCNA作為生物網(wǎng)絡(luò)分析中最具代表性的分析方法[5]，可以多層次多數(shù)量地分析生物系統(tǒng)中基因表達(dá)的信息，已經(jīng)為多物種的基因表達(dá)分析提供了具有臨床意義的結(jié)果。由于目前缺少關(guān)于BLCA的治療效果及預(yù)后的基因?qū)用鏅C(jī)制的系統(tǒng)大數(shù)據(jù)研究。因此，本研究開展了WGCNA方法篩選BLCA預(yù)后及治療效果的相關(guān)基因。

1 材料和方法

1.1 材料

1.1.1 訓(xùn)練數(shù)據(jù)集 本研究從癌癥基因組圖譜(the cancer genome atlas，TCGA)數(shù)據(jù)庫[6]中獲取同時(shí)具有mRNA測序及相應(yīng)臨床數(shù)據(jù)的BLCA患者408例作為訓(xùn)練數(shù)據(jù)集。mRNA測序數(shù)據(jù)經(jīng)過RSEM歸一化和log2轉(zhuǎn)換[7]。選擇年齡、性別、身高、體重、人種、生存時(shí)間、吸煙年數(shù)、治療效果、淋巴結(jié)浸潤、病理分級及TNM分期為臨床表型，以生存時(shí)間和治療效果為主要研究結(jié)果。

1.1.2 驗(yàn)證數(shù)據(jù)集 本研究獲取基因表達(dá)匯編(gene expression omnibus，GEO)數(shù)據(jù)庫中的GSE13507和GSE3167兩個(gè)BLCA數(shù)據(jù)集進(jìn)行驗(yàn)證[8,9]。GSE3167包含9例正常膀胱組織和41例BLCA樣本。GSE13057包含68例正常膀胱或癌旁組織和188例BLCA樣本，并且含有隨訪信息。

1.2 方法

1.2.1 WGCNA 本研究采用R語言WGCNA程序包[10]進(jìn)行WGCNA。一般按照以下步驟進(jìn)行：①離群值的去除：將患者按照基因表達(dá)模式進(jìn)行層次聚類，并將離群樣本剔除；②基因模塊的形成：將基因進(jìn)行K-means聚類，形成模塊；③模塊篩選：將模塊進(jìn)行主成分分析，即各模塊第一主成分與臨床表型進(jìn)行皮爾森相關(guān)分析，得到模塊與表型的相關(guān)系數(shù)即模塊顯著性(module significance，MS)及P值；④模塊內(nèi)重要基因篩選：模塊隸屬度(module membership,MM)>0.5為閾值，以基因?qū)χ委熜Ч闹匾?gene significance,GS)大小排序，篩選模塊中與治療效果顯著相關(guān)的基因。

1.2.2 生存分析 本研究利用R語言Survival程序包進(jìn)行生存分析。將篩選出的基因按表達(dá)量大小繪制生存曲線并進(jìn)行Cox回歸分析[11]，獲得與預(yù)后顯著相關(guān)的基因。

1.2.3 通路活性分析 本研究利用GSCALite(http://bioinfo.life.hust.edu.cn/web/GSCALite/)數(shù)據(jù)庫[12]進(jìn)行腫瘤相關(guān)通路活性的分析。通路包括細(xì)胞凋亡、細(xì)胞周期、DNA損傷修復(fù)、上皮細(xì)胞-間充質(zhì)轉(zhuǎn)化、雄激素受體通路、雌激素受體通路、PI3K/Akt,、Ras/MAPK、RTK及TSC/mTOR。該信息由Rehan Akbani等[13]通過多組學(xué)系統(tǒng)生物學(xué)分析得到。

1.3 統(tǒng)計(jì)學(xué)方法 利用R語言3.4.4統(tǒng)計(jì)軟件。在WGCNA中，相關(guān)分析均采用皮爾森檢驗(yàn)；腫瘤和正常組織的基因表達(dá)量差異使用t檢驗(yàn)；生存分析采用單因素Cox回歸進(jìn)行。P<0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 WGCNA結(jié)果

2.1.1 去除低表達(dá)基因和離群樣本 將在70%及以上樣本中無表達(dá)的基因去除后，408個(gè)樣本進(jìn)行層次聚類，刪去離群的5個(gè)樣本，剩下403個(gè)樣本。樣本聚類及相應(yīng)臨床表型數(shù)據(jù)，圖1A。

2.1.2 基因模塊聚類及模塊篩選 在403個(gè)樣本中進(jìn)行基因模塊聚類，共18 196個(gè)基因聚類為89個(gè)基因模塊(圖1B)。paleturquoise模塊與生存時(shí)間(MS=-0.13，P=0.009)和治療效果顯著相關(guān)(MS=-0.27，P<0.000 1)，但與年齡無關(guān)(MS=-0.058，P>0.05)。因此，本研究選擇paleturquoise模塊做后續(xù)分析。

圖1 WGCNA

2.1.3 模塊內(nèi)重要基因篩選 在paleturquoise模塊中，選取GS值排序前20的基因(表1)，進(jìn)行單因素Cox回歸分析。結(jié)果顯示：8個(gè)基因GPHN、CPSF2、PRMT5、EIF2S1、KLC1、SLC39A9、SRP54及CNIH為預(yù)后相關(guān)因子。

表1 GS值前20基因的信息

2.2 多數(shù)據(jù)集驗(yàn)證

2.2.1 關(guān)鍵基因在腫瘤和正常組織中的表達(dá) 本研究用GEO數(shù)據(jù)庫中的GSE3167和GSE13507分析關(guān)鍵基因在腫瘤和正常組織中的表達(dá)(圖2A)。大多數(shù)基因在腫瘤和正常組織中顯示出了差異表達(dá)。

2.2.2 關(guān)鍵基因與預(yù)后的關(guān)系 將上述8個(gè)關(guān)鍵基因按照表達(dá)量高低分為兩組，并繪制生存曲線。與預(yù)期一致，這8個(gè)基因在TCGA數(shù)據(jù)中均為預(yù)后的不利因素(圖2B)。在GSE13507中趨勢與TCGA一致，除CPSF2和SLC39A9之外，其他基因均為預(yù)后的風(fēng)險(xiǎn)因素(圖2C)。

圖2 關(guān)鍵基因的多數(shù)據(jù)集驗(yàn)證

2.3 關(guān)鍵基因與腫瘤通路活性有關(guān) 通路活性分析顯示GPHN、CPSF2、PRMT5、EIF2S1、KLC1、SLC39A9、SRP54及CNIH在10個(gè)腫瘤相關(guān)通路中均顯示出一定活性(圖3)。幾乎所有的基因都對受體酪氨酸激酶(receptor tyrosine kinase,RTK)有一定活性?；钚暂^強(qiáng)的基因通路關(guān)系還有CPSF2、EIF2S1、PRMT5對細(xì)胞周期，EIF2S1對細(xì)胞凋亡，KLC1對EMT和TSC/mTOR通路。

圖3 通路活性分析

3 討論

本研究對408個(gè)臨床患者樣本信息進(jìn)行WGCNA，構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)，從89個(gè)基因模塊中篩選與患者療效和生存均顯著的paleturquoise模塊，并且該模塊與年齡無關(guān)。通過設(shè)定MM值的閾值，以GS值大小排序，與生存顯著相關(guān)的基因取交集，最終篩選出與預(yù)后及治療效果顯著相關(guān)的GPHN、CPSF2、PRMT5、EIF2S1、KLC1、SLC39A9、SRP54和CNIH基因。通過對上述基因的腫瘤相關(guān)通路活性分析，發(fā)現(xiàn)上述基因都有一定激活RTK通路的作用。GPHN是參與膜蛋白細(xì)胞骨架的微管相關(guān)蛋白[14]；CPSF2主要與mRNA剪接和RNA聚合酶Ⅱ轉(zhuǎn)錄終止有關(guān)[15]；PRMT5編碼屬于甲基轉(zhuǎn)移酶家族的酶，編碼的蛋白質(zhì)催化甲基基團(tuán)轉(zhuǎn)移到氨基酸精氨酸中，包括組蛋白、轉(zhuǎn)錄延伸因子和腫瘤抑制因子p53[16]；EIF2S1是翻譯起始因子EIF2的亞基，與蛋白質(zhì)合成有關(guān)[17]；KLC1編碼驅(qū)動蛋白輕鏈家族的成員，與微管運(yùn)動活性有關(guān)[18]；SLC39A9與金屬離子轉(zhuǎn)運(yùn)及Nrf2途徑有關(guān)[19]；SRP54與mRNA剪接有關(guān)，并且具有GTP酶活性[20]；CNIH是AMPA受體輔助蛋白，涉及TGF-α家族蛋白的選擇性轉(zhuǎn)運(yùn)和成熟[21]。

目前有不少文獻(xiàn)將WGCNA方法運(yùn)用到多種腫瘤的研究中，如腎腫瘤[22]、前列腺癌[23]、肺腺癌[24]、乳腺癌[25]及腸癌[26]等。也有研究關(guān)注膀胱癌的WGCNA，但與本研究的方法及結(jié)果均有所不同。李丹妮[27]的研究關(guān)注女性特異性膀胱癌相關(guān)基因，Deng等[28]的研究關(guān)注腫瘤特異性基因，Li等[29]的研究關(guān)注腫瘤TNM分期。本研究關(guān)注預(yù)后及治療效果，更注重于臨床意義。因?yàn)楸狙芯繉τ陬A(yù)后相關(guān)基因的篩選使用了WGCNA與Cox回歸相結(jié)合的方法，對于模塊的篩選及其他影響因素如年齡進(jìn)行了限制，所以所得結(jié)果更具有穩(wěn)健性。

本研究同樣具有局限性：①只篩選了paleturquoise模塊，可能會漏掉同樣有顯著意義的其他模塊里的基因；②多數(shù)據(jù)集中，這些基因在BLCA和正常組織中的表達(dá)并不完全一致，這可能與樣本的異質(zhì)性以及樣本組織中腫瘤細(xì)胞的純度有關(guān)；③關(guān)于上述基因在BLCA作用的研究還相當(dāng)少，樣本量有限，機(jī)制仍不明確。因此，未來還需要更多的研究來加以驗(yàn)證。

本研究發(fā)現(xiàn)GPHN、CPSF2、PRMT5、EIF2S1、KLC1、SLC39A9、SRP54及CNIH基因與BLCA的療效和預(yù)后均顯著相關(guān)，并且8個(gè)基因都對RTK有一定關(guān)聯(lián)；其中，CPSF2、EIF2S1、PRMT5對細(xì)胞周期有較大相關(guān)性，EIF2S1對細(xì)胞凋亡，KLC1對EMT和TSC/mTOR通路都有較大的關(guān)聯(lián)性。它們有望成為BLCA新的預(yù)后標(biāo)記物或潛在治療靶點(diǎn)，為未來的基礎(chǔ)和臨床研究提供了新的思路。