基于加權基因共表達網絡的膀胱尿路上皮癌預后及療效相關基因篩選

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

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

加權基因共表達網絡分析(weighed gene co-expression network analysis,WGCNA)是一種構建基因共表達網絡的系統生物學算法，該算法以轉錄組數據為分析基礎，在世界生物醫學范圍均得到廣泛的應用[4]。WGCNA作為生物網絡分析中最具代表性的分析方法[5]，可以多層次多數量地分析生物系統中基因表達的信息，已經為多物種的基因表達分析提供了具有臨床意義的結果。由于目前缺少關于BLCA的治療效果及預后的基因層面機制的系統大數據研究。因此，本研究開展了WGCNA方法篩選BLCA預后及治療效果的相關基因。

1 材料和方法

1.1 材料

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

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

1.2 方法

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

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

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

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

2 結果

2.1 WGCNA結果

2.1.1 去除低表達基因和離群樣本 將在70%及以上樣本中無表達的基因去除后，408個樣本進行層次聚類，刪去離群的5個樣本，剩下403個樣本。樣本聚類及相應臨床表型數據，圖1A。

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

圖1 WGCNA

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

表1 GS值前20基因的信息

2.2 多數據集驗證

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

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

圖2 關鍵基因的多數據集驗證

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

圖3 通路活性分析

3 討論

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

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

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

本研究發現GPHN、CPSF2、PRMT5、EIF2S1、KLC1、SLC39A9、SRP54及CNIH基因與BLCA的療效和預后均顯著相關，并且8個基因都對RTK有一定關聯；其中，CPSF2、EIF2S1、PRMT5對細胞周期有較大相關性，EIF2S1對細胞凋亡，KLC1對EMT和TSC/mTOR通路都有較大的關聯性。它們有望成為BLCA新的預后標記物或潛在治療靶點，為未來的基礎和臨床研究提供了新的思路。