基于加權(quán)基因共表達網(wǎng)絡分析篩選結(jié)腸癌關鍵長鏈非編碼RNA及其競爭性內(nèi)源RNA網(wǎng)絡構(gòu)建△

李倩，盧金磊，王會新，崔馨桐，王建喬，侯曉雯，馮旭

沈陽醫(yī)學院公共衛(wèi)生學院，沈陽 110034

結(jié)腸癌是西歐、北美等發(fā)達國家最常見的惡性腫瘤，也是中國最常見的惡性腫瘤之一[1]。由于缺乏早期癥狀，且臨床常用的腫瘤標志物缺乏對早期結(jié)腸癌的診斷效能，大部分患者確診時已處于中晚期，預后較差[2]。近些年研究發(fā)現(xiàn)，長鏈非編碼RNA（long non-coding RNA，lncRNA）在腫瘤的發(fā)生、發(fā)展、預后及轉(zhuǎn)歸中發(fā)揮著重要作用，與結(jié)腸癌發(fā)生發(fā)展相關的lncRNA報道也逐漸增多[3]。加權(quán)基因共表達網(wǎng)絡分析（weighted gene co-expression network analysis，WGCNA）被廣泛用于生物基因研究中，通過聚類的方式更加快捷地找到關鍵基因，同時發(fā)現(xiàn)關鍵基因可能的功能，極大提高了研究速度及準確性[4-6]。本研究通過來自癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數(shù)據(jù)庫和GEO數(shù)據(jù)庫的數(shù)據(jù)，構(gòu)建結(jié)腸癌lncRNA的共表達網(wǎng)絡，篩選得到的lncRNA能夠為進一步研究結(jié)腸癌的潛在發(fā)病機制提供參考，現(xiàn)報道如下。

1 材料與方法

1.1 數(shù)據(jù)下載及處理

從GEO數(shù)據(jù)庫（https://www.ncbi.nlm.nih.gov/geo/）下載GSE126092芯片中的數(shù)據(jù)，包括10對結(jié)腸癌組織及癌旁組織。從TCGA數(shù)據(jù)庫（https//portal.gdc.cancer.gov/）下載結(jié)腸癌轉(zhuǎn)錄組數(shù)據(jù)，其中包括結(jié)腸癌組織482例和正常結(jié)腸組織42例。分析工具主要為R軟件（R×64 4.02版本）及各類R包、Cytoscape（3.8.0版本）以及各類在線數(shù)據(jù)分析網(wǎng)站。

1.2 結(jié)腸癌相關lncRNA的篩選

通過R軟件中的limma程序包，對GEO數(shù)據(jù)進行背景校正、標準化處理以及差異表達分析，篩選標準為：|logFC|≥1.5，校正后P＜0.05。差異分析的結(jié)果用R軟件中pheatmap程序包繪制的火山圖進行可視化分析。

對TCGA數(shù)據(jù)庫下載的結(jié)腸癌組織和正常結(jié)腸組織表達譜數(shù)據(jù)進行WGCNA分析，首先進行離群值的篩選，隨后進行軟閾值的確定，使用R軟件自帶的層次聚類函數(shù)hclust進行聚類分析，使用不同的顏色標記聚類分析中的模塊。模塊與樣本信息進行相關性分析，從中選擇與結(jié)腸癌相關性最高的模塊，獲取該模塊基因進行后續(xù)分析。

對GEO數(shù)據(jù)中的差異表達lncRNA和WGCNA性狀相關模塊中的lncRNA取交集，獲取關鍵lncRNA，進行后續(xù)分析。

1.3 競爭性內(nèi)源RNA（competing endogenous RNA，ceRNA）調(diào)控網(wǎng)絡的構(gòu)建

對上述獲得的lncRNA進行ceRNA網(wǎng)絡的構(gòu)建，使用Starbase（http://starbase.sysu.edu.cn/）預測lncRNA的靶向miRNA，使用miRDB數(shù)據(jù)庫（http://mirdb.org/）和Targetscan數(shù)據(jù)庫（http://www.targetscan.org/）預測miRNA的靶基因mRNA?；谏鲜龊Y選出的lncRNA、miRNA、mRNA，采用Cytoscape（3.8.0版本）構(gòu)建并繪制ceRNA調(diào)控網(wǎng)絡。

1.4 構(gòu)建蛋白質(zhì)-蛋白質(zhì)相互作用（protein-protein interaction，PPI）網(wǎng)絡

使用String數(shù)據(jù)庫構(gòu)建PPI網(wǎng)絡。輸入基因集為 mRNA；種屬選擇為 Homo sapiens；combined score≥0.7。使用Cytoscape（3.8.0版本）軟件可視化PPI數(shù)據(jù)。

1.5 基因功能分析

使用DAVID在線數(shù)據(jù)庫（https://david.ncifcrf.gov/）進行mRNA的基因本位（Gene Ontology，GO）功能分析和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）富集分析。

2 結(jié)果

2.1 結(jié)腸癌相關lncRNA的篩選

2.1.1 GEO 中差異表達lncRNA的篩選GEO數(shù)據(jù)庫GSE126092芯片中共篩選出322個差異表達的lncRNA，包含113個上調(diào)基因和209個下調(diào)基因。（圖1）

圖1 GEO數(shù)據(jù)庫GSE126092芯片中差異表達的lncRNA火山圖

2.1.2 WGCNA 分析結(jié)果 經(jīng)樣本聚類分析后刪除15個離群樣本。為使得鄰接函數(shù)滿足無尺度網(wǎng)絡的條件，選取β=3進行后續(xù)分析，此時共表達網(wǎng)絡接近為無尺度網(wǎng)絡。根據(jù)β=3進行切割設置得到基因聚類樹，每個模塊最少基因數(shù)目設置為30，得到7個lncRNA模塊（圖2A）。對模塊與樣本特征進行相關性分析，最終確定綠色模塊（cor=0.85，P＜0.05）為與結(jié)腸癌相關性最高的模塊（圖2B）。對GEO中差異表達的lncRNA和綠色模塊中的lncRNA取交集，最終獲得6個結(jié)腸癌的關鍵lncRNA，分別為鋅指NFX1結(jié)構(gòu)1反義RNA1（zinc finger NFX1-type containing 1 antisense RNA 1，ZFAS1），β1,3-半乳糖基轉(zhuǎn)移酶5反義RNA1（beta-1,3-galactosyltransferase5antisenseRNA 1，B3GALT5-AS1），細胞色素P450家族1亞家族B成員1反義RNA1（cytochrome P450 family 1 subfamily B member 1 antisense RNA 1，CYP1B1-AS1），二肽基肽酶樣10反義RNA1（dipeptidyl peptidase like 10 antisense RNA 1，DPP10-AS1），VPS9包含域1反義RNA1（VPS9 domain containing 1 antisense RNA 1，VPS9D1-AS1）和細胞周期蛋白依賴性激酶抑制因子2B反義RNA1（cyclin dependent kinase inhibitor 2B antisense RNA 1，CDKN2B-AS1）。

圖2 WGCNA 分析結(jié)果

2.2 構(gòu)建ceRNA調(diào)控網(wǎng)絡

預測出與6個關鍵lncRNA可能相互作用的24個miRNA，以及24個miRNA可能的靶基因mRNA共346個，構(gòu)建了lncRNA介導的ceRNA網(wǎng)絡圖。（圖3）

圖3 構(gòu)建結(jié)腸癌的lncRNA-miRNA-mRNAceRNA網(wǎng)絡圖

2.3 PPI網(wǎng)絡構(gòu)建

構(gòu)建PPI網(wǎng)絡圖鑒定ceRNA網(wǎng)絡中mRNA蛋白質(zhì)間的相互作用關系，發(fā)現(xiàn)一些聯(lián)合評分比較高的mRNA，分別為：雌激素受體1（estrogen receptor 1，ESR1）、小窩蛋白1（caveolin1，CAV1）、間質(zhì)-上皮細胞轉(zhuǎn)化因子（mesenchymal-epithelial transition factor，MET）、鈣黏蛋白相關蛋白β1（cadherinassociated protein beta 1，CTNNB1）、磷脂酰肌醇轉(zhuǎn)移蛋白 3（phosphatidylinositol transfer protein 3，PITPNM3）和趨化因子 18（chemokine ligand 18，CCL18）。（圖4）

圖4 PPI網(wǎng)絡圖

2.4 基因功能分析

對346個mRNA進行GO功能分析和KEGG富集分析。GO功能分析結(jié)果顯示，生物功能主要集中在DNA模板轉(zhuǎn)錄調(diào)控、RNA聚合酶Ⅱ基因啟動子的轉(zhuǎn)錄調(diào)控和RNA聚合酶Ⅱ啟動子轉(zhuǎn)錄負向調(diào)控等；細胞功能主要集中在細胞核、突觸、神經(jīng)細胞體、突觸后密集區(qū)和微管相關復合體；分子功能主要集中在核酸結(jié)合、金屬離子結(jié)合、DNA結(jié)合等。KEGG富集分析結(jié)果顯示，基因主要富集在癌癥蛋白聚糖、磷脂酰肌醇-3-羥激酶（phosphatidylinositol 3-hydroxy kinase，PI3K）-蛋白激酶 B（protein kinase B，PKB，又稱AKT）信號通路、Rap1信號通路和局部粘連等。（圖5、圖6）

圖5 GO功能分析

圖6 KEGG富集分析

3 討論

結(jié)腸癌發(fā)生與社會環(huán)境、高脂肪飲食、遺傳等密切相關，具有發(fā)病率高、轉(zhuǎn)移率高、治愈率低等特點[7]。因此，非常有必要在分子水平上開發(fā)新的生物標志物和潛在靶點以預防和治療結(jié)腸癌。WGCNA可以通過系統(tǒng)繪制個體生物網(wǎng)絡互作圖精準找出與研究相關的核心基因，極大提高了研究速度及準確性[4,8]。因此，本研究通過構(gòu)建結(jié)腸癌WGCNA共表達網(wǎng)絡，尋找與結(jié)腸癌具有密切關聯(lián)性的lncRNA。

本研究從GEO數(shù)據(jù)庫中共篩選出322個差異基因，從TCGA數(shù)據(jù)庫中共獲得1688個lncRNA的表達矩陣，進行WGCNA構(gòu)建后，綠色模塊為與結(jié)腸癌相關性最高的模塊，對GEO數(shù)據(jù)中差異表達的lncRNA和TCGA綠色模塊中的lncRNA取交集后最終獲得6個結(jié)腸癌的關鍵lncRNA。構(gòu)建的ceRNA網(wǎng)絡提示其在結(jié)腸癌中的可能作用機制，但仍需進一步實驗驗證。

已有研究表明這6種lncRNA在腫瘤的發(fā)生發(fā)展中發(fā)揮重要作用。ZFAS1定位于人類染色體20q13，研究發(fā)現(xiàn)ZFAS1與結(jié)腸癌的分化程度、T分期及N分期有關，高表達ZFAS1是結(jié)腸癌預后不良的危險因素[9-10]。馮偉[11]發(fā)現(xiàn)，胃癌患者血清B3GALT5-AS1表達上調(diào)，可能作為潛在的胃癌輔助診斷及預后監(jiān)測的生物標志物。另有研究發(fā)現(xiàn)，DPP10-AS1、VPS9D1-AS1和CDKN2B-AS1均具有促進肺癌細胞增殖的作用，可促進肺癌惡性進展[12-14]。雖然現(xiàn)有研究表明B3GALT5-AS1、CYP1B1-AS1、DPP10-AS1、VPS9D1-AS1和 CDKN2B-AS1與腫瘤的發(fā)生發(fā)展有關，但并未有研究表明它們與結(jié)腸癌相關，因此如將其作為一種診斷指標，仍需進一步研究以提供更可靠的依據(jù)。

綜上所述，本研究利用GEO數(shù)據(jù)庫和TCGA數(shù)據(jù)庫以及WGCNA方法篩選出與結(jié)腸癌可能相關的6個lncRNA，分別為ZFAS1、B3GALT5-AS1、CYP1B1-AS1、DPP10-AS1、VPS9D1-AS1和 CDKN2B-AS1，并且構(gòu)建了相關的ceRNA調(diào)控網(wǎng)絡，為進一步探索結(jié)腸癌的機制研究提供了依據(jù)。