基于轉錄組測序構建食管腺癌預后相關的內源競爭RNA網絡

陳超 李冬 童國俊 胡四平

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

1 對象和方法

1.1 對象 選擇截至2019年7月從TCGA獲取78例食管腺癌組織和13例正常食管組織的lncRNA、miRNA、mRNA表達量數據，以及食管腺癌組織患者性別、年齡、病理分期等臨床數據。食管腺癌患者男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的差異倍數（fold change，FC）以2為底的對數的絕對值（|log2FC|）均為≥1；（2）mRNA 的|log2FC|≥0.6；（3）P＜0.05。以lncRNA、miRNA、mRNA的表達量的中位值將食管腺癌樣本分為各基因所對應的高表達組和低表達組，并統計每個差異表達的lncRNA、miRNA、mRNA的基因數。

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

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

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）結果。對富集分析結果根據標準化富集分數（normalized enrichment score，NES）進行排序，對各組中NES排名前5的KEGG通路結果進行可視化。

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

2 結果

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

2.2 生存分析與權重基因共表達網絡分析結果 （1）在所有表達有差異的lncRNA、miRNA、mRNA中，生存分析結果提示與預后相關的RNA數目分別有81、20、374個。（2）由這些與預后相關的 lncRNA、miRNA、mRNA組成的權重基因共表達網絡中，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又存在共表達關系的預后相關的mRNA有7個，分別為CDK1、TPX2、細絲蛋白A（filamin A，FLNA），肌肉相關受體酪氨酸激酶（muscle associated receptor tyrosine kinase，MUSK）、蛋白激酶camp依賴的Ⅱ型調節亞基β（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 食管腺癌預后相關的ceRNA網絡構建

2.3.1 預后相關的ceRNA網絡 見圖1。

由圖1可見，在該ceRNA網絡中AC131009.3、AC087388.1、miR-29c-3p、CDK1、TPX2 互相之間存在內源競爭網絡關系，即構成了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調控網絡。

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

2.3.2 ceRNA網絡中各基因表達水平與生存時間的關系 見圖2。

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

圖2 ceRNA網絡中各基因表達水平與生存時間的關系

2.3.3 ceRNA網絡中各基因在食管腺癌組織、正常食管組織中的表達差異 見圖3。

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

圖3 ceRNA網絡中各基因在食管腺癌組織、正常食管組織中的表達差異

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

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

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

2.4 GSEA中的可能與腫瘤相關的KEGG富集結果CDK1、TPX2基因參與的可能與腫瘤預后相關的KEGG見表1。

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

表1 CDK1、TPX2基因參與的可能與腫瘤預后相關的KEGG

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

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

3 討論

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

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

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

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