肺腺癌中綜合分析競爭性內(nèi)源RNA網(wǎng)絡預測潛在的預后lncRNA

黃銀龍，連超群，孫康，張志強，鄧嬌嬌，張靜*

（1.蚌埠醫(yī)學院生命科學學院，安徽 蚌埠 233030；2.檢驗醫(yī)學院）

當前肺腺癌是促使癌癥患者死亡的因素之一［1］，盡管肺腺癌的治療方法不斷的改進，但是肺腺癌的預后效果還是相對較差［2］，肺腺癌晚期患者的癌細胞會發(fā)生侵襲和轉(zhuǎn)移，導致對患者的治療效果較差［3］。肺腺癌的發(fā)生和發(fā)展機制難于探索，預防和治療肺腺癌還是人類面臨的巨大挑戰(zhàn)之一。因此更好地理解肺腺癌的侵襲和轉(zhuǎn)移的基礎機制，并且鑒定出與肺腺癌相關的新診斷方法和預后生物靶標是目前的主要任務。

構(gòu)建lncRNA、miRNA和mRNA的競爭性內(nèi)源RNA（ceRNA）網(wǎng)絡。lncRNA和mRNA可以直接與miRNA結(jié)合，因而lncRNA可以作為內(nèi)源競爭性分子與miRNA結(jié)合從而間接地調(diào)控mRNA的表達［4］。這種間接的lncRNA-mRNA調(diào)節(jié)關系的作用已被證明與多種癌癥發(fā)生發(fā)展直接相關。例如，已證明DGCR5通過競爭性結(jié)合hsa-mir-22-3p促進Bcl-2的上調(diào)，從而促進肺腺癌細胞的增殖和轉(zhuǎn)移［5］。但是缺乏可供查詢的與腫瘤相關的lncRNAmiRNA-mRNA關系的ceRNA網(wǎng)絡數(shù)據(jù)庫，進而難以理解許多癌癥病理學基礎。

癌癥基因組圖譜（TCGA）是一個大型公共數(shù)據(jù)庫，存儲30多種人類癌癥相對應的數(shù)據(jù)，包含患者的大量臨床病理信息［6］。鑒于其強大的特性，TCGA非常適合數(shù)據(jù)挖掘工作，探索腫瘤發(fā)展的基礎機制。該數(shù)據(jù)庫已被用于肝癌［7］、乳腺癌［8］、胃癌［9］和舌鱗狀細胞癌［10］等腫瘤類型ceRNA網(wǎng)絡的構(gòu)建。此類ceRNA網(wǎng)絡主要用于強調(diào)復雜遺傳關系和鑒定潛在的預后或治療性生物標志物。

本研究主要比較肺腺癌患者組織樣品中的lncRNA、miRNA和mRNA與正常組織樣品之間的差異表達，并關注這些差異表達基因與肺腺癌分期的關系。然后運用肺腺癌不同階段之間通用的差異RNAs（lncRNA、miRNA、mRNA）構(gòu)建與肺腺 癌 相 關 的ceRNA網(wǎng) 絡。通 過miRcode、TargetScan、miRDB和miRTarBase數(shù)據(jù)庫預測了不同差異RNAs之間的關系，以預測它們之間的相互作用關系。獲得由170個lncRNA、14個miRNA和66個mRNA組成的ceRNA網(wǎng)絡。由此產(chǎn)生的網(wǎng)絡關系可以提供有價值的并且能快速識別肺腺癌相關的差異基因。此外，我們通過Cox回歸模型等分析方法生成基于6個lncRNA的風險評分模型，該模型能夠預測肺腺癌患者的生存結(jié)果，從而突出了該研究方法的價值。

1 材料與方法

1.1 數(shù)據(jù)與預處理從TCGA數(shù)據(jù)庫中下載肺腺癌患者的lncRNA、miRNA、mRNA及其相關的臨床數(shù)據(jù)［11-12］。一共獲得594個樣本，其中包括535個癌癥樣本和59個正常樣本。運用perl［13］提取樣本信息并轉(zhuǎn)化成矩陣文件，使用Ensemble數(shù)據(jù)庫把RNA矩陣數(shù)據(jù)中的基因序列號轉(zhuǎn)換成基因名［14］。對包含3種RNA類型、生存時間、生存狀態(tài)和分期的樣本進行分析，確保每個樣本中包含miRNA、lncRNA、mRNA數(shù)據(jù)以及用于預后分析的臨床數(shù)據(jù)。

1.2 篩選差異基因根據(jù)臨床數(shù)據(jù)的腫瘤分期人為設定為前期（Ⅰ～Ⅱ期）、后期（Ⅲ～Ⅳ期），提取前期、后期、正常樣本中的RNAs數(shù)據(jù)。應用R軟件中的edgeR包［15-16］將前期樣本和后期樣本分別與正常樣本進行RNAs差異分析，設定P＜0.01，|logFC|＞1作為差異分析的標準［17］。將得到lncRNA、miRNA和mRNA的差異基因用R軟件包gplots繪制火山圖［18］，使用R軟件中的VennDiagram包獲得通用差異RNAs［19］。

1.3 ceRNA網(wǎng)絡的構(gòu)建運用miRcode［20］預測差異lncRNA在miRNA中的目標基因，使用miRDB［21］、miRTarBase［22］、TargetScan［23］預 測miRNA在mRNA的靶基因，從而獲得差異lncRNAmiRNA-mRNA的網(wǎng)絡關系，應用Cytoscape軟件對該網(wǎng)絡關系進行可視化［24］。

1.4 mRNA生物學功能與通路分析基因本體（gene ontolontolgy，GO）和京都基因和基因組百科全書（kyoto encyclopedia of genes and genomes，KEGG）對ceRNA網(wǎng)絡中的mRNA進行富集分析，從而研究肺腺癌發(fā)生的分子基礎。利用R軟件中的clusterProfiler包［25］對mRNA進行GO分析［26］（生物過程、細胞組成、分子功能）和KEGG分析［27］（生物學通路），設定篩選標準為P＜0.05。

1.5 生存分析鑒于lncRNA在生命調(diào)節(jié)中獨特的生物學意義，因此對ceRNA網(wǎng)絡中的差異lncRNA進行了單變量Cox回歸模型分析［28］，進而揭示其在肺癌患者生存時間之間的關聯(lián)。應用R軟件中的survival包對其分析，計算危險比（HR）和相應的95%可信區(qū)間（95%CI），篩選標準為P＜0.05。

1.6 構(gòu)建風險模型使用R軟件中的survminer包對單因素回歸分析得到的顯著性lncRNA進行風險模型構(gòu)建，然后計算每個樣品的風險值，并根據(jù)生存中位值對每例患者進行分組，分為高風險組與低風險組。應用R軟件對分組情況進行可視化得到風險圖和狀態(tài)圖。通過R軟件對每個樣本進行生存曲線分析，最后運用ROC曲線確定我們生存曲線模型的準確性。

1.7 獨立預后分析通過單變量Cox分析，初步評定了特定的臨床信息（年齡、性別、腫瘤分級、TN分期、病理分期）與肺腺癌患者的關系。通過單變量Cox分析，將P＜0.05的變量納入多變量分析中作為候選變量。在多變量Cox回歸模型分析中［29］，評定HR和相應的95%CI，篩選標準為P＜0.05。

1.8 PPI網(wǎng)絡構(gòu)建和預測評估通過回歸分析得到的lncRNA構(gòu)建與肺腺癌患者預后相關的子網(wǎng)絡，使用STRING對該網(wǎng)絡中的mRNA進行蛋白質(zhì)-蛋白質(zhì)互作分析。隨后應用Cytoscape的插件CytoHubba從lncRNA-mRNA網(wǎng)絡中挑選出9個關鍵基因。利用R軟件對6個lncRNA和9個mRNA相關性分析。運用在線軟件Kaplan-Meier對顯著性相關的核心基因進行生存分析，評定標準為P＜0.01，|R|＞0.3。利用在線GEPIA軟件對8對lncRNA-mRNA關系對中的基因在癌癥和正常組織中的表達情況進行分析。

2 結(jié)果

2.1 lncRNA、mRNA、miRNA差異表達鑒定為了識別肺腺癌分期相關的差異RNAs，將下載的樣本分為正常樣本（n=59），肺腺癌早期樣本（Ⅰ～Ⅱ期，n=403）和肺腺癌晚期樣本（Ⅲ～Ⅳ期，n=111）。然后將早期肺腺癌樣本和晚期肺腺癌樣本分別與正常樣本比較，獲得差異lncRNA、miRNA和mRNA，篩選條件為|logFC|＞1和P＜0.01（圖1）。比較了早期和晚期樣本的差異lncRNA、miRNA和mRNA，得 到 通 用 的4 847個mRNA、2 070個lncRNA和232個miRNA。這些差異表達的重疊基因可能同時導致癌癥的發(fā)生與發(fā)展，因此為了更好地理解這些已識別的差異RNAs（lncRNA、miRNA、mRNA）之間的調(diào)控關系，需要更進一步對lncRNA-miRNA和miRNA-mRNA的關系對進行研究。

圖1 肺腺癌相關的差異RNAs的鑒定

2.2 鑒定lncRNA-miRNA相互作用關系利用重疊差異RNAs鑒定lncRNA-miRNA之間的關系。詳細的流程見圖2。首先使用miRcode數(shù)據(jù)庫查找2 070個lncRNA的目標基因miRNA。確定了87個能代表差異miRNA的miRNA。最后預測到190個lncRNA和36個miRNA之間存在相互作用關系。

2.3 鑒定miRNA-mRNA相互作用關系使用重疊差異miRNA鑒定miRNA-mRNA之間的相互作用關系。詳 細 流 程 見 圖2。首 先 利 用miRDB、miRTarBase和TargetScan數(shù) 據(jù) 庫 查 找232個 差 異miRNA的靶基因mRNA，然后將得到的靶基因（mRNA）與差異mRNA比較取交集。最后我們預測到66個mRNA和14個miRNA之間存在聯(lián)系。

圖2 ceRNA網(wǎng)絡構(gòu)建的流程圖

2.4 ceRNA網(wǎng)絡構(gòu)建利用lncRNA-miRNA和miRNA-mRNA關系對，最終得到190個lncRNA和66個mRNA。又通過lncRNA-miRNA-mRNA關系對刪除部分lncRNA-miRNA與miRNA-mRNA關系對。最終得到170個lncRNA、14個miRNA和66個mRNA，設定閾值為P＜0.05。運用Cytoscape軟件構(gòu)建了ceRNA網(wǎng)絡（圖3）。

圖3 構(gòu)建與肺腺癌相關的ceRNA網(wǎng)絡

2.5 富集分析通過GO和KEGG數(shù)據(jù)庫對ceRNA網(wǎng)絡中的66個mRNA潛在的生物學功能和生物學通路進行富集分析。發(fā)現(xiàn)mRNA在24種生物學過程中顯著富集（adjustedP＜0.05）。mRNA主要參與RNA聚合酶Ⅱ特異性、DNA結(jié)合轉(zhuǎn)錄激活劑活性、L-谷氨酸跨膜轉(zhuǎn)運蛋白活性、Wnt信號體等生物學功能（圖4A）。得到mRNA在4種KEGG通路顯著性富集（adjustedP＜0.05），包括AGE-RAGE信號通路、轉(zhuǎn)錄失調(diào)、氨基酸的生物合成（圖4B）。這些結(jié)果表明AGE-RAGE信號通路、轉(zhuǎn)錄失調(diào)通路對肺腺癌發(fā)展起到重要的調(diào)節(jié)作用。

圖4 ceRNA網(wǎng)絡中66個差異表達mRNA的功能富集分析

2.6 與肺腺癌生存相關的lncRNA為了確定與肺腺癌生存顯著性相關的lncRNA，對ceRNA網(wǎng)絡中的lncRNA進行了單因素回歸模型分析，最后得到7個與肺腺癌患者的整體生存期顯著性相關的lncRNA，DGCR5（HR=1.000 398，95%CI：1.000 01～1.000 787，P＜0.05），LINC00536（HR=1.021 126，95%CI：1.004 467～1.038 062，P＜0.012 739），ATG9B（HR=0.999 15，95%CI：0.998 46～1.000 787，P＜0.015 752），F(xiàn)GF14-IT1（HR=0.868 353，95%CI：0.766 564～0.983 659，P＜0.026 489），ANO1-AS2（HR=1.006 849，95%CI：1.002 042～1.011 679，P＜0.005 18），ST7-AS2（HR=1.001 523，95%CI：1.000 077～1.002 972，P＜0.038 976），SYNPR-AS1（HR=0.994 116，95%CI：0.988 899～0.999 362，P＜0.027 961）

2.7 構(gòu)建lncRNA相關風險評分系統(tǒng)鑒于lncRNA明確的表達方式和在ceRNA網(wǎng)絡中的頂層調(diào)節(jié)因子作用，因此我們認為其可以作為肺腺癌診斷和預后治療的理想的生物標志物。使用R軟件survminer包對7個顯著性lncRNA構(gòu)建Cox回歸模型，得到6個lncRNA（DGCR5、LINC00536、ATG9B、FGF14-IT1、ANO1-AS2、SYNPR-AS1）（P＜0.05）與風險評分模型，風險值：0.000 553 593×DGCR5＋0.019835865×LINC00536＋（-0.000941518）×ATG9B＋（-0.109 168 446）×FGF14-IT1＋0.007 943 741×ANO1-AS2＋（-0.005 158 638）×SYNPR-AS1。使用上述公式確定了所有樣本的患者的風險得分，并根據(jù)風險得分是否大于1.245 336（風險中位值）將患者分成高風險組（n=169）和低風險組（n=169）。運用R軟件對分組情況進行可視化得到風險圖（圖5A）和狀態(tài)圖（圖5B），可以清晰地得出患者樣本在高低風險組的分布狀況。根據(jù)分組情況，利用R軟件繪制高低風險組的生存曲線（圖5C），可以得出風險模型把患者分為高低風險組，高風險組的死亡率高于低風險組。并利用R軟件中的survival包繪制ROC曲線（圖5D）驗證生存曲線的準確性，證明生存曲線具有一定的可靠性。運用單因素獨立預后分析，驗證肺腺癌患者中與總生存時間相關的臨床因素（圖5E），又運用多變量Cox回歸分析，探索與總生存時間相關的臨床獨立因素（圖5F），單因素分析發(fā)現(xiàn)腫瘤分期（Stage）、淋巴結(jié)（N）、遠處轉(zhuǎn)移（T）和風險模型值可以作為獨立預后因子，多因素分析發(fā)現(xiàn)風險模型值可以作為獨立預后因子。

圖5 lncRNA相關風險評分系統(tǒng)

2.8 PPI網(wǎng)絡構(gòu)建和預后評估根據(jù)風險模型得到的6個lncRNA構(gòu)建了預后子網(wǎng)絡（圖6A），并使用該網(wǎng)絡中的mRNA構(gòu)建PPI網(wǎng)絡（圖6B）。從該網(wǎng)絡中 挑 選 了9個 核 心 基 因（FGF2、CCNB1、TGFBR2、ALDOA、BDNF、COL1A1、MCM4、PBK、SLC1A1）用于進一步研究。根據(jù)lncRNAmiRNA-mRNA的關系，理論上我們認為lncRNA可以調(diào)節(jié)mRNA的表達。因此對6個lncRNA和9個mRNA進行了相關性分析，結(jié)果得到8對lncRNAmRNA（|R|＞0.30），具 體 是FGF2、CCNB1、ALDOA、COL1A1和MCM4與6個lncRNA顯著相關（圖7）。然后利用Kaplan-Meier數(shù)據(jù)庫繪制FGF2、CCNB1、ALDOA、COL1A1和MCM4的總生存時間曲線（圖8）。利用GEPIA數(shù)據(jù)庫分析了8對lncRNA-mRNA關系對中的基因（FGF2、CCNB1、ALDOA、COL1A1、MCM4、DGCR5、ATG9B）在癌癥和正常組織中的表達情況（圖9）。最終可以得出基因FGF2、CCNB1、ALDOA、COL1A1、MCM4、DGCR5、ATG9B可以作為肺腺癌治療的獨立預后因子。

圖6 ceRNA網(wǎng)絡和PPI網(wǎng)絡的構(gòu)建

圖7 對6個與風險評分相關的lncRNA和9個核心基因進行相關性分析

圖8 基因表達水平對肺腺癌患者總生存時間的影響

圖9 lncRNA-mRNA關系對中基因在肺腺癌和正常樣本的表達情況

3 討論

近年來，肺癌是患病率和死亡率明顯增加的一種惡性腫瘤［30］，非小細胞肺癌又是肺癌的主要類型之一，占所有病例約85%，而肺腺癌在非小細胞肺癌中的占比最高，約40%［31］。雖然手術切除可以一定程度上改善早期肺癌患者，但是其它原因會導致肺癌復發(fā)或轉(zhuǎn)移。因此進一步了解肺癌的發(fā)生和發(fā)展機制是發(fā)現(xiàn)新的治療手段的重要環(huán)節(jié)。隨著高通量測序技術的發(fā)展，揭示了lncRNA在眾多生物學過程中發(fā)揮著重要的調(diào)控作用。越來越多的證據(jù)表明，lncRNA在ceRNA網(wǎng)絡中發(fā)揮重要作用，并且ceRNA網(wǎng)絡中的靶基因?qū)τ诜伟┗颊叩念A后及耐藥性的治療可能具有重要作用［32-34］。本研究發(fā)現(xiàn)6個lncRNA中有一個是以前沒有報道過的（FGF14-IT1），有5個是以前已經(jīng)報道 過 的（DGCR5、LINC00536、ATG9B、ANO1-AS2、SYNPR-AS1）。其中DGCR5、SYNPR-AS1是已報道并與肺癌相關的lncRNA，而LINC00536、ATG9B、ANO1-AS2是與肺癌不相關，但與其它癌癥緊密相關的lncRNA。DGCR5是一種人們比較了解的致癌因子之一，它通過競爭性結(jié)合miR-218-5p促進肺癌的遷移與侵襲［35］。這種lncRNA已被證明可以參與其他類型的腫瘤細胞的形成、增殖和轉(zhuǎn)移，包括宮頸癌［36］、胰腺癌［37］、膽囊癌［38］和肝癌［39-40］，在這些腫瘤類型中，它與患者較差的預后密切相關。LINC00536可以作為膀胱癌［41］、乳腺癌［42-43］和卵巢癌［44］等癌癥的生物標志物以及在其中具有預后作用。ATG9B被證實在腎細胞癌［45］、肝細胞癌［46］、乳腺癌［47］的增殖和遷移的發(fā)展中發(fā)揮重要作用。SYNPR-AS1也被驗證在肺癌的發(fā)展中發(fā)揮重要的作用［48］。因此，我們的ceRNA網(wǎng)絡既能夠成功識別先前驗證過的肺腺癌相關的lncRNA，如DGCR5、SYNPR-AS1，又能夠識別先前報道極少的FGF14-IT1。進一步分析還發(fā)現(xiàn)與6個lncRNA相關的核心基因（CCNB1、ALDOA、COL1A1、MCM4）也參與了腫瘤的發(fā)生與發(fā)展［49-52］。

綜上所述，我們利用不同疾病分期的肺腺癌患者的組織樣本數(shù)據(jù)生成一個包含lncRNAmiRNA-mRNA關系對的ceRNA網(wǎng)絡。這個全面的網(wǎng)絡不僅有利于對肺腺癌的基礎機制的研究，也突出了對未來診斷和治療的潛在靶點研究，強調(diào)了將lncRNA作為肺腺癌患者預后和治療潛在生物標志物的價值。由于以往構(gòu)建肺腺癌相關lncRNA網(wǎng)絡的研究較少，因此我們的研究存在一定的局限性。值得注意的是，我們通過分析4個核心基因的存活率，間接地對6個lncRNA生存率進行了外部驗證。此外在風險評分模型中，6個lncRNA的機制作用還沒有被直接驗證。因此，還需要進一步的工作來驗證和擴大我們的發(fā)現(xiàn)。