基于轉(zhuǎn)錄組測序構(gòu)建食管腺癌預(yù)后相關(guān)的內(nèi)源競爭RNA網(wǎng)絡(luò)

陳超 李冬 童國俊 胡四平

食管腺癌已成為發(fā)病率最高的食管惡性腫瘤，也是生長速度最快、最具侵襲性的惡性腫瘤之一[1-2]。由于目前食管腺癌的預(yù)后仍極不理想，探索食管腺癌的發(fā)生、發(fā)展及預(yù)后影響因素是食管腺癌治療的基礎(chǔ)。微小非編碼 RNA（microRNA，miRNA）可通過識別 mRNA上的miRNA應(yīng)答元件抑制mRNA的翻譯或?qū)е耺RNA降解，同時，長鏈非編碼RNA（long non-coding RNA，lncRNA）上也存在miRNA應(yīng)答元件并競爭性結(jié)合miRNA[3]。因此，lncRNA、miRNA、mRNA之間形成的內(nèi)源競爭RNA（competing endogenous RNA，ceRNA）網(wǎng)絡(luò)是食管癌預(yù)后機制研究的重要突破點，可能為食管腺癌患者的治療提供新的方法。癌癥和腫瘤基因圖譜計劃（The Cancer Genome Atlas，TCGA）數(shù)據(jù)庫中包含了大量的腫瘤標本的數(shù)據(jù)，包括數(shù)百例食管惡性腫瘤患者的轉(zhuǎn)錄組測序數(shù)據(jù)及相對應(yīng)的臨床基本信息。然而，關(guān)于食管腺癌中的ceRNA網(wǎng)絡(luò)及其影響預(yù)后的機制目前仍知之甚少。本研究提取食管腺癌組織的轉(zhuǎn)錄組數(shù)據(jù)，利用生物信息分析手段構(gòu)建食管腺癌預(yù)后相關(guān)的ceRNA網(wǎng)絡(luò)，以進一步探索食管腺癌預(yù)后影響機制，現(xiàn)將結(jié)果報道如下。

1 對象和方法

1.1 對象 選擇截至2019年7月從TCGA獲取78例食管腺癌組織和13例正常食管組織的lncRNA、miRNA、mRNA表達量數(shù)據(jù)，以及食管腺癌組織患者性別、年齡、病理分期等臨床數(shù)據(jù)。食管腺癌患者男67例，女11例，年齡 68.4（58.0～77.1歲），生存期為 2.15（1.01～4.88）年；按腫瘤病理分期，處于Ⅰ期、Ⅳ期各10例（12.8%），Ⅱ期25例（32.1%），Ⅲ期 33例（42.3%），13例正常食管組織TCGA未提供臨床信息。

1.2 方法

1.2.1 差異基因篩選 采用R軟件分析食管腺癌組織和正常食管組織中存在差異表達的lncRNA、miRNA、mRNA。差異基因的篩選條件：（1）lncRNA、miRNA的差異倍數(shù)（fold change，F(xiàn)C）以2為底的對數(shù)的絕對值（|log2FC|）均為≥1；（2）mRNA 的|log2FC|≥0.6；（3）P＜0.05。以lncRNA、miRNA、mRNA的表達量的中位值將食管腺癌樣本分為各基因所對應(yīng)的高表達組和低表達組，并統(tǒng)計每個差異表達的lncRNA、miRNA、mRNA的基因數(shù)。

1.2.2 生存分析及權(quán)重基因共表達網(wǎng)絡(luò)分析 用Kaplan-Meier法分析食管腺癌組織中各lncRNA、miRNA、mRNA的表達量的多少與生存時間的關(guān)系，若P＜0.05即說明該基因為食管腺癌預(yù)后相關(guān)基因。用R軟件對差異表達且與預(yù)后相關(guān)的lncRNA和miRNA、miRNA和mRNA構(gòu)建無尺度網(wǎng)絡(luò)，以無尺度網(wǎng)絡(luò)指數(shù)＞0.85為符合無尺度網(wǎng)絡(luò)條件并確定軟閾值。計算并獲得lncRNA與miRNA、miRNA與mRNA的權(quán)重基因共表達網(wǎng)絡(luò)。

1.2.3 預(yù)后相關(guān)的ceRNA網(wǎng)絡(luò)構(gòu)建及可視化 以miRNA為中心，在權(quán)重基因共表達網(wǎng)絡(luò)中篩選與miRNA存在共表達關(guān)系的mRNA，同時逆向篩選可吸附miRNA并與預(yù)后相關(guān)的lncRNA，構(gòu)建lncRNA/miRNA/mRNA內(nèi)源競爭網(wǎng)絡(luò)，并在Cytoscape軟件中可視化；再次采用生存分析評估ceRNA網(wǎng)絡(luò)中各基因的表達水平與預(yù)后的關(guān)系，分析ceRNA網(wǎng)絡(luò)中各種RNA在食管腺癌組織、正常食管組織中的表達差異，以及miRNA與lncRNA、mRNA表達的相關(guān)性。

1.2.4 基因集富集分析（gene set enrichment analysis,GSEA） 將食管腺癌組織分別分為細胞周期蛋白激酶1（cyclin dependent kinase 1，CDK1）高表達組、CDK1 低表達組、微管成核因子（microtubule nucleation factor,TPX2）高表達組、TPX2低表達組，隨后進行GSEA研究京都基因與基因組百科全書（Kyoto encyclopedia of genes and genomes，KEGG）結(jié)果。對富集分析結(jié)果根據(jù)標準化富集分數(shù)（normalized enrichment score，NES）進行排序，對各組中NES排名前5的KEGG通路結(jié)果進行可視化。

1.3 統(tǒng)計學處理 數(shù)據(jù)統(tǒng)計分析及可視化采用R軟件（3.6.1版）、Cytoscape軟件（3.7.2版）。計量資料如呈正態(tài)分布以表示，組間比較采用t檢驗，否則以四分位數(shù)表示，組間比較采用非參數(shù)檢驗。分類變量采用例數(shù)（百分率）表示，生存資料采用log-rank檢驗，相關(guān)性分析采用Pearson相關(guān)，P＜0.05為差異有統(tǒng)計學意義。

2 結(jié)果

2.1 差異基因篩選結(jié)果 與正常食管組織比較，食管腺癌組織中符合篩選條件被分到高表達組的lncRNA、miRNA、mRNA分別有678、131、1 862個，被分到低表達組的 lncRNA、miRNA、mRNA分別有 774、59、2 825個。

2.2 生存分析與權(quán)重基因共表達網(wǎng)絡(luò)分析結(jié)果 （1）在所有表達有差異的lncRNA、miRNA、mRNA中，生存分析結(jié)果提示與預(yù)后相關(guān)的RNA數(shù)目分別有81、20、374個。（2）由這些與預(yù)后相關(guān)的 lncRNA、miRNA、mRNA組成的權(quán)重基因共表達網(wǎng)絡(luò)中，lncRNA共有7個，分別為 AC131009.3、AC024337.2、AC245407.1、BDNF-AS、AL136172.1、AC087388.1、AC010261.2；miRNA 共有 3個，分別為 miR-125b-5p、miR-143-3p、miR-29c-3p，而與這3個miRNA又存在共表達關(guān)系的預(yù)后相關(guān)的mRNA有7個，分別為CDK1、TPX2、細絲蛋白A（filamin A，F(xiàn)LNA），肌肉相關(guān)受體酪氨酸激酶（muscle associated receptor tyrosine kinase，MUSK）、蛋白激酶camp依賴的Ⅱ型調(diào)節(jié)亞基β（protein kinase cAMP-dependent type Ⅱ regulatory subunit beta，PRKAR2B）、非受體5型蛋白酪氨酸磷酸酶（protein tyrosine phosphatase，non-receptor type 5，PTPN5）、血管活性腸肽受體 2（vasoactive intestinal peptide receptor 2，VIPR2）。

2.3 食管腺癌預(yù)后相關(guān)的ceRNA網(wǎng)絡(luò)構(gòu)建

2.3.1 預(yù)后相關(guān)的ceRNA網(wǎng)絡(luò) 見圖1。

由圖1可見，在該ceRNA網(wǎng)絡(luò)中AC131009.3、AC087388.1、miR-29c-3p、CDK1、TPX2 互相之間存在內(nèi)源競爭網(wǎng)絡(luò)關(guān)系，即構(gòu)成了ACO87388.1/miR-29c-3p/CDK1、ACO87388.1/miR-29c-3p/TPX2、AC131009.3/miR-29c-3p/CDK1、AC131009.3/miR-29c-3p/TPX2 的ceRNA調(diào)控網(wǎng)絡(luò)。

圖1 食管腺癌的預(yù)后ceRNA網(wǎng)絡(luò)圖，菱形為lncRNA，三角形為miRNA，橢圓形為mRNA；深灰色為與正常食管組織比較，在食管腺癌組織中高表達，淺灰色為食管腺癌組織中低表達（TPX2為微管成核因子，CDK1為細胞周期蛋白激酶1，F(xiàn)LNA為細絲蛋白A，MUSK為肌肉相關(guān)受體酪氨酸激酶，PRKAR2B為蛋白激酶camp依賴的Ⅱ型調(diào)節(jié)亞基β，PTPN5為非受體5型蛋白酪氨酸磷酸酶，VIPR2為血管活性腸肽受體2）

2.3.2 ceRNA網(wǎng)絡(luò)中各基因表達水平與生存時間的關(guān)系 見圖2。

由圖 2 可見，AC131009.3、AC087388.1、CDK1、TPX2在食管腺癌組織中的表達水平越高，患者的生存時間越短，生存率越低；miR-29c-3p在食管腺癌組織中的表達水平越高，患者的生存時間越長，生存率越高。即miR-29c-3p在食管腺癌中可能是腫瘤抑制基因，AC131009.3、AC087388.1、CDK1、TPX2 在食管腺癌中可能是促腫瘤生長基因。

圖2 ceRNA網(wǎng)絡(luò)中各基因表達水平與生存時間的關(guān)系

2.3.3 ceRNA網(wǎng)絡(luò)中各基因在食管腺癌組織、正常食管組織中的表達差異 見圖3。

由圖3可見，與正常食管組織比較，在食管腺癌組織中 ACO87388.1、AC131009.3、CDK1、TPX2 的表達水平更高，而miR-29c-3p在食管腺癌組織中的表達水平更低，差異有統(tǒng)計學意義。ceRNA網(wǎng)絡(luò)中各基因在食管腺癌組織、正常食管組織中的表達差異趨勢與食管腺癌患者生存分析結(jié)果相符。

圖3 ceRNA網(wǎng)絡(luò)中各基因在食管腺癌組織、正常食管組織中的表達差異

2.3.4 miR-29c-3p 表達量與 ACO87388.1、AC131009.3、CDK1、TPX2表達量的相關(guān)性分析 見圖4。

由圖4可見，在共表達分析中，miR-29c-3p的表達量與ceRNA網(wǎng)絡(luò)中的lncRNA（AC087388.1、AC131009.3）、mRNA（CDK1、TPX2）的表達量存在明顯的Pearson相關(guān)，相關(guān)系數(shù)分別為-0.473、-0.564、-0.655、-0.647。

圖4 miR-29c-3p表達量與ACO87388.1、AC131009.3、CDK1、TPX2表達量的相關(guān)性分析

2.4 GSEA中的可能與腫瘤相關(guān)的KEGG富集結(jié)果CDK1、TPX2基因參與的可能與腫瘤預(yù)后相關(guān)的KEGG見表1。

由表1可見，KEGG顯示的是CDK1高表達組、CDK1低表達組、TPX2高表達組、TPX2低表達組中NSE排名前5的KEGG通路；CDK1基因可能參與腫瘤的免疫、細胞周期等生物學過程有關(guān)，TPX2基因主要與腫瘤的代謝等。對表1進行可視化，CDK1和TPX2基因參與的KEGG結(jié)果見圖5（插頁）。

表1 CDK1、TPX2基因參與的可能與腫瘤預(yù)后相關(guān)的KEGG

圖5 細胞周期蛋白激酶1（CDK1）和微管成核因子（TPX2）基因參與的KEGG結(jié)果（a：CDK1；b：TPX2）

由圖5可見，GSEA結(jié)果提示 CDK1、TPX2基因在食管腺癌中可參與多種高級生物學功能，以及各生物學功能在CDK1高表達組、CDK1低表達組、TPX2高表達組、TPX2低表達組中的分布情況。

3 討論

食管腺癌是常見的食管惡性腫瘤，腫瘤細胞對常規(guī)抗癌治療的耐藥率很高，晚期食管癌的5年生存率僅為0.9%[4]。食管腺癌表觀遺傳現(xiàn)象異常引起的生物學功能改變是導(dǎo)致腫瘤患者預(yù)后不佳的重要因素之一，而基因表達水平的改變是表現(xiàn)遺傳現(xiàn)象建立的核心。基因表達水平的異常是導(dǎo)致腫瘤預(yù)后差異的起源，分子間的信號傳導(dǎo)是預(yù)后改變的重要機制。Zhang等[5]就曾利用TCGA中的食管腺癌數(shù)據(jù)并采用權(quán)重基因共表達網(wǎng)絡(luò)分析方法探索食管腺癌中的預(yù)后基因。近年來已有大量文獻表明，具有差異表達的lncRNA、miRNA可能在食管癌患者的預(yù)后中扮演者重要的作用[6-10]。非編碼RNA的生物學作用主要是通過復(fù)雜的互作網(wǎng)絡(luò)最終作用于具有蛋白編碼功能的mRNA而發(fā)揮作用[11]。因此，研究食管腺癌中l(wèi)ncRNA、miRNA、mRNA之間的競爭網(wǎng)絡(luò)來探索腫瘤的生物學價值可能為腫瘤的治療提供新的治療手段[12]。

生物信息分析技術(shù)是一種高效率、精準確率、高可信度的數(shù)據(jù)挖掘技術(shù)，已被《Nature》等知名期刊認可[13-14]。本研究中筆者通過生物信息分析篩選食管腺癌中存在差異表達且與預(yù)后相關(guān)的lncRNA、miRNA、mRNA，再挖掘到以miRNA為中心的ceRNA網(wǎng)絡(luò)。ceRNA網(wǎng)絡(luò)構(gòu)建方式是多樣化的，例如可以通過靶基因預(yù)測工具或權(quán)重基因共表達網(wǎng)絡(luò)分析等方法來實現(xiàn)[15-17]，本研究中的ceRNA是通過分析相應(yīng)樣本中l(wèi)ncRNA、miRNA、mRNA的共表達關(guān)系的權(quán)重系數(shù)獲得的，研究結(jié)果可能更貼近實際，可靠性更強。該ceRNA網(wǎng)絡(luò)的核心是miR-29c-3p，多項研究表明miR-29c-3p在ceRNA中可以抑制卵巢癌生長[18]，在室管膜瘤中可以改變腫瘤細胞對化療藥物的敏感性等[19]。在食管腺癌中，Craig等[20]學者報道m(xù)iR-29c-3p出現(xiàn)了異常表達，與本研究結(jié)果一致。本研究共表達分析結(jié)果提示miR-29c-3p與ACO87388.1、AC131009.3、CDK1、TPX2 呈負相關(guān)關(guān)系，與CDK1、TPX2的相關(guān)性最為密切（Cor＞-0.6）。然而，食管腺癌中關(guān)于miR-29c-3p在ceRNA網(wǎng)絡(luò)中的作用及其機制報道甚少。在現(xiàn)有研究中，ceRNA網(wǎng)絡(luò)中的ACO87388.1、AC131009.3與惡性腫瘤的相關(guān)性尚不明確，但是已有證據(jù)表明CDK1、TPX2是重要的腫瘤預(yù)后相關(guān)分子[21-23]。因此，進一步地研究ceRNA網(wǎng)絡(luò)在食管腺癌中的生物學作用顯得極為迫切。

lncRNA、miRNA最終是通過影響mRNA的翻譯或使mRNA降解而發(fā)揮生物學作用，因此CDK1、TPX2是本研究的核心分子。基于CDK1、TPX2基因RNA表達量進行的GSEA結(jié)果認為，本ceRNA網(wǎng)絡(luò)可能參與了食管腺癌的多種生物學功能，例如影響食管腺癌的離子通到、細胞代謝、腫瘤微環(huán)境等途徑從而調(diào)控腫瘤細胞的增殖、侵襲等。在惡性腫瘤中，離子通道的轉(zhuǎn)運因改變了腫瘤的代謝而影響了患者的預(yù)后，是惡性腫瘤靶向治療的靶點[24-25]。此外，腫瘤微環(huán)境可通過表達異常的分子建立強大的腫瘤免疫逃逸系統(tǒng)，在食管鱗癌中腫瘤的微環(huán)境可影響腫瘤的生長[26]，在食管腺癌也可能存在類似的機制。

綜上所述，食管腺癌是一種高度惡性的腫瘤，本研究發(fā)現(xiàn) ACO87388.1、AC131009.3、miR-29c-3p、CDK1、TPX2是影響食管腺癌患者的編碼、非編碼基因，以miR-29c-3p為中心的ceRNA可能參與了食管腺癌腫瘤細胞的發(fā)生、發(fā)展過程，在未來食管腺癌的預(yù)后分析中存在重要的研究價值。本研究僅局限于生物信息的分析研究，研究結(jié)論尚需進一步的基礎(chǔ)實驗、臨床研究去證實，但研究結(jié)果為未來的研究探索提供了思路。