基于加權基因共表達網絡挖掘卵巢癌相關基因△

李康梅，陳泯锜，戴明明，陳秀榕，黎明星，何國珍#

廣西中醫藥大學1基礎醫學院，2賽恩斯新醫藥學院，3研究生院，南寧 530200

卵巢癌是婦女生殖系統致死率最高的惡性腫瘤，發病率和病死率均較高，并呈逐年上升趨勢[1-2]。數據表明，中國每年有14 080人死于卵巢癌[3]。盡管腫瘤的治療不斷發展，但仍有2/3以上的卵巢癌患者未能早期發現和確診，生存率仍無明顯改善，且其發病機制和起源尚未清楚[4]。因此，研究卵巢癌的發生發展機制，尋找相關的樞紐基因對卵巢癌的早期診治及預后判斷具有一定意義。

加權基因共表達網絡分析（weighted gene coexpression network analysis，WGCNA）是一種基于高通量基因表達譜的算法，廣泛應用于各種疾病的基因共表達網絡識別中，以揭示基因間的相關性，尋找顯著相關的基因模塊。借助WGCNA算法構建基因模塊，可挖掘基因模塊信息，對相關基因模塊中的樞紐基因進行研究。本研究基于WGCNA對美國國立生物技術信息中心（National Center for Biotechnology Information，NCBI）的基因表達綜合數據庫（gene expression omnibus，GEO）下載的卵巢癌基因表達譜數據集GSE18520進行挖掘相關基因，構建卵巢癌基因共表達網絡，識別卵巢癌相關基因模塊及模塊內樞紐（Hub）基因，為進一步探索卵巢癌的發生發展機制提供一種新思路。

1 材料與方法

1.1 數據下載和數據預處理

從GEO數據庫（https://www.ncbi.nlm.nih.gov/geo/）下載基因表達譜數據集GSE18520，共53個卵巢癌樣本和10個正常樣本，以篩選卵巢癌的差異表達基因（differentially expressed gene，DEG）。對下載的數據進行預處理，包括提取樣本信息、構建基因表達矩陣以及將探針名轉化為基因名等。

1.2 對卵巢癌進行DEG篩選

使用“limma”R包篩選卵巢癌樣本和正常樣本之間的DEG。將調整后的P＜0.05和|log2 fold change（FC）|＞1作為篩選條件。

1.3 卵巢癌中DEG的PPI網絡

使用 STRING（http://www.string-db.org/）數據庫評估基因列表背后的相互關系。同時，對網絡進行MCODE分析，以degree cut-off=5，node score cut-off=0.2，k-core=2，maximum depth=100為參考值選擇PPI網絡的核心模塊。

1.4 卵巢癌中DEG的共表達分析

安裝R軟件的WGCNA包，為DEG構建一個無標度的共表達網絡，選取合適的加權系數β，將相關矩陣轉化為鄰接矩陣。建樹的方法采用動態混合剪切法，將相異度作為距離測度，設定最小模塊尺寸為30，進行模塊識別并繪制基因樹狀圖。

1.5 核心模塊及核心基因的識別

WGCNA中的每一個模塊都在STRING數據庫中進行分析，通過模塊內統計分析，根據每個基因與該模塊的相關系數即模塊隸屬度（module membership，MM）≥0.8且定義基因的顯著性（gene significance，GS）P＜0.05篩選出核心基因。

1.6 功能注釋和通路分析

為了研究大量基因的功能注釋，使用“Clusterpro-filer”軟件包對樞紐模塊進行了基因本體（gene ontology，GO）功能富集分析和京都基因以及基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）通路富集分析。

1.7 樞紐基因的驗證和生存分析

利用基因表達譜交互分析數據庫（gene expression profiles interactive analysis database，GEPIA，http://GEPIA.cancer-pku.cn/index.html）對核心基因進行進一步驗證和生存分析。以P＜0.01在GEPIA數據庫中繪制樞紐基因的生存曲線。

1.8 統計學方法

采用R語言平臺Rstuio對卵巢癌臨床特征和基因過表達模塊相關系數進行計算。

2 結果

2.1 卵巢癌的DEG篩選

在R軟件的“limma”軟件包中，對53個卵巢癌樣本和10個正常樣本的GSE18520基因表達譜進行分析。共檢測到2474個DEG（上調864個，下調1610個）。（圖1）

圖1 GSE18520基因表達譜中所有DEG的火山圖

2.2 樞紐模塊的鑒定

根據在MCODE中評估的分數得出的前3個簇作為核心模塊。其中MCODE 1有40個節點和592條邊，在這些集群中得分最高（圖2）。

2.3 篩選軟閾值

共表達網絡符合無尺度網絡，即出現連接度為k的節點的對數lgk與該節點出現的概率的對數lg[P（k）]呈負相關，且相關系數應＞0.8。使用R軟件WGCNA包構建權重共表達網絡，使用分析包自動選擇的軟閾值計算得到軟閾值β=20（scale freeR^2=0.89）。

2.4 WGCNA 的構建和關鍵模塊的識別

使用WGCNA R包，可以通過平均連鎖聚類將具有高度相關表達模式的DEG放入模塊中，并挖掘出8個模塊（圖3A）。設置閾值為0.25以合并相似模塊（圖3B），并相應生成6個模塊。無法納入任何模塊的基因放入灰色模塊，在后續分析中剔除。其他5個模塊分別用藍色、棕色、黑色、粉色和藍綠色表示（圖3A）。GS和MM分析在藍色和藍綠色模塊中顯示出非常顯著的相關性（圖3C）。

圖2 在MOCDE中評估的PPI網絡中得分最高的樞紐模塊

2.5 樞紐模塊的GO功能富集分析和KEGG通路富集分析

對藍色和藍綠色樞紐模塊進行GO功能富集分析以及KEGG通路富集分析。GO功能富集分析分別有生物過程（biological process，BP）涉及的調節神經遞質水平和細胞組成（cellular component，CC）涉及的細胞皮層、核外膜、細胞器外膜、有機外膜（表1）。KEGG通路富集分析主要有癌癥蛋白聚糖、視黃醇的新陳代謝、黏著斑、WNT信號通路、肌動蛋白骨架的調節、代謝途徑（表2）。

2.6 樞紐基因的驗證和生存分析

圖3 WGCNA的構建和關鍵模塊的識別

表1 藍色和藍綠色樞紐模塊GO功能富集分析

表2 藍色和藍綠色樞紐模塊KEGG富集通路分析

通過GEPIA數據庫對藍色模塊和藍綠色模塊的482個核心基因進行驗證和生存分析。根據模塊相關性大于0.9進行篩選，得到以下8個基因：AC104667.3、BTD、DDX26B、K+通道子族B成員1（K+channel subfamily B member 1，KCNB1）、前列腺素F2α受體（prostaglandin F2 receptor，PTGFR）、腦糖原磷酸化酶（brain glycogen phosphorylase，PYGB）、runt相關轉錄因子1（runt-related transcriptionfactor1，RUNX1T1）、TMEM91。其中AC104667.3、PYGB在卵巢癌中表達上調，BTD、DDX26B、KCNB1、PTGFR、RUNX1T1、TMEM91在卵巢癌中表達下調。AC104667.3、DDX26B的高表達與更好的總生存率有關，BTD、KCNB1、PTGFR、PYGB、RUNX1T1、TMEM91的高表達與預后不良有關（圖4）。

3 討論

本研究主要通過生物信息學分析GSE18520（53個卵巢癌樣本和10個正常樣本）的基因表達譜中卵巢癌發生及發展過程中的樞紐基因。研究結果顯示，2474個基因在腫瘤組織中有差異表達。在這2474個基因中，有864個上調基因和1610個下調基因。利用Cytoscape軟件的MCODE插件選擇PPI網絡中最重要的前3個模塊作為樞紐模塊，并對最重要的樞紐模塊（MCODE1）進行GO功能富集分析以及KEGG通路富集分析的進一步研究。

圖4 不同基因表達水平的卵巢癌患者的生存曲線

為了在卵巢癌中發現具有潛在價值的基因，在藍色模塊和藍綠色模塊的482個樞紐基因中根據相關系數大于0.9篩選出8個樞紐基因（AC104667.3、BTD、DDX26B、KCNB1、PTGFR、PYGB、RUNX1T1、TMEM91）。Tan 等[5]研究表明，BTD在宮頸癌中表達上調，作為宮頸癌的潛在基因，其意義還需進一步研究。DDX26B是一種RNA解旋酶，是不同的RNA進行代謝所必需的[6]。劉志榮等[7]研究表明DDX26B基因在早期大腸癌和腺瘤中不表達，而在正常大腸組織中表達。KCNB1通道的氧化是一種神經脆弱的機制，這種機制在脊椎動物的大腦中普遍存在[8]。如小鼠衰老過程中KCNB1通道的中度氧化可以通過影響突觸功能影響神經元網絡[9]。Yu等[10]通過使用藥物直接撞擊KCNB1通道氧化激活的毒性路徑，改善小鼠顱腦損傷（traumatic brain injury，TBI）的結果，從而闡明了脊椎動物神經毒性的機制，為人類創傷性腦損害提供一種新的治療辦法。研究認為PTGFR和神經生長因子（nerve growth factor，NGF）蛋白在子宮腺肌病痛經的發生和發展中起作用，血清NGF蛋白水平及血漿PTGFR蛋白水平與子宮腺肌病痛經程度呈正相關，通過抑制其表達可能為子宮腺肌病痛經提供新的治療線索[11-12]。研究表明，PTGFR激活可能誘發子宮內膜修復移植PTGS-2基因表達和刺激血管內皮生長因子（vascular endothelial growth factor，VEGF）、結締組織生長因子（connective tissue growth factor，CTGF）、轉化生長因子β1（transforming growth factor β1，TGF-β1），引發蛋白激酶C（protein kinase C，PKC）信號通路激活和基因表達[13]。PYGB目前已被報道參與多種腫瘤的發生，糖原代謝與許多其他關鍵代謝途徑交叉，在腫瘤微環境中發揮著關鍵作用，抑制糖原代謝可能是一種有前途的抗癌治療策略[14-15]。Zhou等[16]研究表明，PYGB在卵巢癌組織中表達上調，高水平的PYGB表達與卵巢癌患者預后不良顯著相關，這與本研究結果一致。另外敲除PYGB后可顯著抑制卵巢癌細胞的增殖、侵襲和遷移，并發現miRNA-133a-3p/PYGB軸在卵巢癌治療中是一個重要的診斷指標和治療靶點，可能成為未來治療卵巢癌的靶標[16]。Sun等[4]研究發現，EPB41L4A-AS2通過與miRNA-103a聯合誘導RUNX1T1表達可導致卵巢癌細胞的遷移和侵襲，證明了EPB41L4AAS2、miRNA-103a和RUNX1T1的相互作用在卵巢癌的發生和發展中起關鍵作用，但關于miRNA-103a和RUNX1T1的表達如何影響卵巢癌的發生和發展還需要進一步的研究。Yeh等[17]研究也表明了RUNX1T1在卵巢癌細胞系中表達下調，可以抑制卵巢癌細胞的生長，通過卵巢癌中異常的TGF-β/SMAD家族成員4（SMAD family member 4，SMAD4）信號建立表觀遺傳沉默，但RUNX1T1在卵巢癌中的作用機制尚不明確，還需進一步研究。另外，RUNX1T1蛋白也作為診斷肝轉移的一個有前途的生物標志物[4]。而AC104667.3、TMEM91目前尚未有相關研究。

綜上所述，本研究通過運用一系列合理的生物信息學手段，結合文獻報道以及本研究結果，發現8個基因均與卵巢癌的發生和發展相關。而PYGB和RUNX1T1基因與卵巢癌的生存及預后密切相關，可能通過細胞周期相關通路影響卵巢癌的發生、進展及預后，為卵巢癌的進一步研究提供依據及參考。