通過加權基因共表達網絡分析(WGCNA)研究食管腺癌的關鍵基因

姚 歌，王新帥，張晨曦，李新陽，王子明，王鈺鯤，張衛國

食管癌是發生在食管上皮組織的惡性腫瘤，是全世界范圍內癌癥相關死亡的主要原因之一，占所有惡性腫瘤的2%，死亡率高達90%，發生率和死亡率在惡性腫瘤中分別居第六位和第四位[1]。食管癌主要有兩種組織學亞型：鱗狀細胞癌和腺癌，都是高度侵襲性腫瘤。最新的流行病學數據表明，全球79%的食管鱗癌發生在東南部和中亞，而46%的食管腺癌發生于北歐、西歐、北美洲和大洋洲[2]。在過去的幾年中，發達國家腺癌發病率有所上升，而鱗癌的發病率下降。同時食管癌仍然是世界范圍內最難治愈的惡性腫瘤之一，由于診斷時多為晚期，只有不到50%的食管癌患者有手術切除的機會。近年來隨著手術、化療和放療技術的發展，食管癌的臨床治療有了很大進展，但5 a總生存率仍然很低。分子靶向治療是在基因分子水平上，針對已經明確的致癌位點來設計相應的治療藥物，藥物進入體內會特異地結合致癌位點發生作用，使腫瘤細胞特異性壞死，較少對正常組織影響。隨著腫瘤分子靶向治療在肺癌、結直腸癌等疾病的成功應用，加上靶向藥物本身選擇性高、抗腫瘤活性強優點，分子靶向治療已成為人們試圖解決食管癌治療問題的新途徑。但目前食管癌可應用的靶點較少，因此尋找新的分子靶點顯得尤為重要。加權基因共表達網絡分析(weighted gene co-expression network analysis,WGCNA)是鑒定表型性狀的基因模塊和關鍵基因的理想方法，已有研究采用此方法對乳腺癌、肺癌、結直腸癌等惡性腫瘤進行分析，但尚未有應用WGCNA對食管腺癌進行分析的研究。因此，本研究采用WGCNA算法識別與食管腺癌發展高度相關的基因模塊，并從中篩選出樞紐基因，為進一步發現新的生物標志及食管腺癌診斷和治療的潛在靶點提供重要的理論依據。現報道如下。

1 材料與方法

1.1 數據來源及數據預處理基因表達譜數據(GSE92396)及其臨床信息來源于美國國家生物信息中心下的GEO數據庫，平臺編號為GPL6244，共納入22例樣本的基因表達數據用于分析，包含9例正常食管組織樣本、12例食管腺癌組織及1例介于兩者之間的組織樣本基因表達譜數據。在進行WGCNA前，對下載的矩陣信息進行了預處理：利用平臺信息將探針名轉化為基因名，將對應多個探針的基因表達量作log2處理及標準化處理后作為該基因的最終表達量。最終獲得行名為樣本名，列名為基因名的矩陣用于后續共表達網絡的構建。

1.2 數據篩選在R語言下運行WGCNA包，進行共表達網絡的構建，為獲取有效的共表達網絡，計算各基因在所有樣本中表達量的方差，取方差前25%的基因用于共表達網絡的構建。對樣本進行聚類分析，檢測并去除離群異常值。

1.3 構建基因共表達網絡通過選擇一個合適的加權系數β(軟閾值)使構建網絡中的基因之間的連接服從無尺度網絡分布，并利用基因之間的相關系數構建分層聚類樹，聚類樹的不同分支代表不同的基因模塊，不同顏色代表不同的模塊。然后基于基因的加權相關系數，將基因按照表達模式進行分類，將模式相似的基因歸為一個模塊，從而將基因通過基因表達模式歸類成不同的模塊來進行下一步分析。最后應用此系數將相關矩陣轉換為鄰接矩陣，進一步轉化為拓撲重疊矩陣(Topological Overlap Matrix,TOM)，并隨機選取400個基因做TOM熱圖來證明模塊之間的高度獨立性以及每個模塊中基因表達的相對獨立性。

1.4 模塊與臨床特征關聯分析通過WGCNA算法計算各模塊基因和樣本構成的矩陣(Module Eigengene,ME)與臨床性狀的皮爾森相關系數及其P值，通過皮爾森相關系數衡量不同模塊與臨床性狀之間的關系，并取相關系數最高的模塊用于后續分析。然后利用基因顯著性(Gene Significance,GS)和模塊顯著性(Module Significance,MS)計算與樣本種類相關模塊的表達模式。GS為每個基因的表達量與某一臨床信息的皮爾森相關系數，MS為模塊中所有基因的GS取絕對值均值。為了篩選目的模塊中的樞紐基因，計算目的模塊中每一基因與該模塊ME的相關系數，即模塊隸屬度(module membership,MM)，用以衡量目的模塊內的某一基因從屬于該模塊的程度。

1.5 樞紐基因的尋找和分析通過對各個模塊和樣本臨床信息進行關聯分析后做出相應熱圖，選取相關性最高的模塊，將該模塊中的基因信息導入到Cytoscape(3.7.1版本)軟件中，使用Cytoscape將模塊中的基因信息可視化為網絡，并通過連接度對模塊中基因進行排序，篩選出前7個基因用于進一步分析。最后把篩選出的基因用GEPIA數據庫做生存分析，觀察其對生存預后的影響。

1.6 統計學處理本研究使用R統計軟件(3.6版)和WGCNA包進行統計計算。本文采用的食管組織的相關臨床特征和各個共表達模塊的ME間相關系數的計算是基于R語言平臺Rstuio下實現。WGCNA即加權基因共表達網絡分析，可被應用于尋找功能相似的基因。WGCNA使用了一種軟閾值(soft thresholding)定義一個權重值，判斷每一組基因連接的可能性。在此概念下，形成加權共表達網絡。WGCNA使用一個以生物學意義作為基準的評判標準——無尺度拓撲學準則。

2 結果

2.1 數據來源及數據預處理從GEO數據庫中得到行名為22個樣本，列名為33 297個核酸探針的矩陣。利用平臺文件GPL6244將探針名轉換為基因名，并對其基因表達量作log2處理及標準化處理，最終得到包含個基因的樣本基因表達譜矩陣，作為共表達網絡構建的輸入文件。

2.2 數據篩選本研究選取方差前25%的14 124個基因用于共表達網絡的構建。通過構建22個食管組織樣本的14 124個基因的層次聚類樹(圖1)，發現無明顯離群樣本。

圖1 22個食管組織樣本系統聚類圖

2.3 構建基因共表達網絡模塊更符合無尺度特征，選取β值為4用于構建共表達網絡(圖2)。通過動態剪切樹法進行模塊識別，最終得到45個模塊(圖3)。隨機選取400個基因做TOM熱圖(圖4)，描述基因之間的拓撲重疊矩陣，確定每個模塊獨立存在。

圖2 加權基因共表達網絡分析軟閾值(β)的確定

圖3 基因聚類樹狀圖

圖4 拓撲重疊矩陣(TOM)圖

2.4 模塊與臨床特征關聯分析對各個模塊和樣本臨床信息進行關聯分析，從模塊和性狀熱圖中可以發現深紅色(darkred)模塊與食管腺癌相關程度最高(圖5)。計算深紅色模塊中基因GS和MM相關系數(cor=0.77，P=5.6e-25)進一步驗證此結果的可信度(圖6)。結果深紅色模塊和食管腺癌的關聯程度最大，因此后續研究選取深紅色模塊篩選樞紐基因。

圖5 臨床特征與模塊特征相關性熱圖

圖6 深紅色模塊基因模塊隸屬度與基因顯著性相關的散點圖

圖7 矩陣網絡熱圖

2.5 尋找并分析樞紐基因進一步分析發現深紅色模塊共有121個基因，通過WGCNA對這 121個基因進行連通性分析后導入Cytoscape軟件做進一步分析。使用Cytoscape軟件對目標模塊的網絡進行可視化(圖8)，并依據連接度選出樞紐基因(圖9)，結果顯示深紅色模塊中的樞紐基因為VASP、BEX4、KDELR3、POLR3B、MYO15B、FBXO27、ABCC3。最后將樞紐基因導入GEPIA數據庫中做生存分析(圖10)。

圖8 深紅色模塊中基因網絡可視化

3 討論

食管癌是世界上最為常見的消化系統惡性腫瘤之一，中國食管癌的發病率和病死率居世界第一位[3]。食管癌惡性程度較高，且早期癥狀不典型，導致早期診斷困難，臨床就診患者多為中晚期，預后極差。隨著食管外科手術技術的提高和放射治療、化學治療技術的不斷創新，食管癌可切除率以及患者生存率有了一定提高，但化學治療藥物毒副作用較大，很多患者無法耐受，且對化學治療藥物有耐藥趨勢[4]。因此，不斷有研究嘗試食管癌生物靶向治療及免疫治療手段，探索相關分子作用的靶點。傳統的基因共表達模式的研究通常通過計算任意兩個基因之間的相關系數，設立硬閾值來衡量兩個基因是否具有相似表達模式，缺乏說服力。WGCNA通過在高通量數據中識別出功能相關或表達相似的基因構成的模塊，從生物功能整體考慮基因功能及其聯系，彌補了傳統方法的缺陷。通過使用這種方法，研究者不僅可以發現某模塊內基因間的聯系，還可以關注該模塊與其他模塊內基因的相關性。此外，通過將臨床信息與模塊相關聯，還可進一步獲得與臨床特征相關的基因，有助于建立疾病某臨床特征的研究基礎[5]。

本研究采用WGCNA對GEO數據庫中的22個食管組織樣本基因表達譜進行分析，最終獲得了45個共表達模塊。將共表達模塊與臨床信息相關聯，從臨床特征與模塊特征相關性熱圖中得到深紅色模塊中的基因與食管腺癌組織相關性系數最高(0.66，P=7e-04)，且與食管正常組織相關性系數最低(-0.59，P=0.004)，深紅色模塊中的基因與食管腺癌的形成高度相關。將深紅色模塊中121個基因間的連通性信息導入到Cytoscape軟件中，進一步識別出7個可能在食管腺癌中發揮關鍵作用的樞紐基因，分別為VASP、POLR3B、BEX4、KDELR3、ABCC3、FBXO27、MYO15B。VASP作為一種肌動蛋白結合蛋白，參與一系列依賴細胞骨架重塑和細胞極性的過程，如軸突導向、片長軸突和絲狀軸突動力學、血小板激活和細胞遷移，VASP與絲狀肌動蛋白的形成有關，可能在細胞黏附和運動中起著廣泛的作用[6]；ABCC3是MDR相關蛋白家族中的一員,能轉運谷胱甘肽、硫酸鹽、膽汁酸和葡糖苷酸聚合物等多種內源性和外源性化合物，被發現在多個惡性腫瘤中有高表達[7]；BEX4通過解除微管動力學和染色體完整性的調控，可抵抗凋亡和細胞死亡，導致非整倍體的獲得，同時增加了腫瘤的增殖潛能和生長[8]；POLR3B在感應和限制細胞內細菌和DNA病毒感染方面起著關鍵作用，在先天免疫反應中起核質和胞質DNA傳感器的作用[9]；KDELR3是編碼KDEL內質網蛋白保留受體家族成員，在酵母和動物細胞中，通過從順式高爾基體或高爾基前腔室中不斷獲取內質網腔內的可溶性蛋白，可以實現內質網腔內可溶性蛋白的保持[10]；FBXO27為SCF (SKP1-CUL1-F-box蛋白)型E3泛素連接酶復合物的底物識別組分，能夠識別和結合變性糖蛋白[11]。MYO15B在食管癌及其他腫瘤中的作用及功能尚未見報道，其作用需要進一步驗證。

圖9 樞紐基因網絡可視化

P<0.05，差異有統計學意義。圖10 利用GEPIA數據庫做樞紐基因的生存分析

將7個樞紐基因導入到GEPIA數據庫中做生存分析，發現VASP、ABCC3、MYO15B基因高表達患者總體生存率低于低表達患者，且P<0.05，說明VASP、ABCC3、MYO15B可能為食管腺癌的致癌基因。BEX4基因高表達患者總體生存率高于低表達患者，BEX4可能為食管腺癌的抑癌基因。因此這些基因可能是影響食管腺癌生存和預后的關鍵基因，可能作為食管腺癌診斷和治療的潛在靶點，對進一步理解食管腺癌的分子機制可起到一定幫助。需要進一步對這些基因的作用進行單獨的分析和實驗驗證。