葡萄膜惡性黑色素瘤轉(zhuǎn)移相關(guān)的非編碼RNA表達(dá)譜及競爭性 內(nèi)源RNA調(diào)控網(wǎng)絡(luò)分析

陳曉云，楊偉敏，鄧小茜，季嫻，肖偉

(中山大學(xué)中山眼科中心，眼科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，廣東省眼科視覺科學(xué)重點(diǎn)實(shí)驗(yàn)室，廣州 510060)

葡萄膜惡性黑色素瘤(uveal melanoma，UM)是成人最常見的眼內(nèi)原發(fā)性惡性腫瘤，好發(fā)于高加索人和西班牙裔。初診治療后，約一半U(xiǎn)M患者會最終發(fā)生遠(yuǎn)處轉(zhuǎn)移[1]。由于缺乏有效的治療手段，大多數(shù)遠(yuǎn)處轉(zhuǎn)移的患者在發(fā)現(xiàn)轉(zhuǎn)移灶后的存活時(shí)間通常少于一年。目前，學(xué)者們已經(jīng)建立了幾種腫瘤基因表達(dá)譜、臨床特征和/或細(xì)胞遺傳學(xué)標(biāo)志物預(yù)測轉(zhuǎn)移風(fēng)險(xiǎn)的模型[2]。最近，癌癥基因組圖譜(The Cancer Genome Atlas，TCGA)項(xiàng)目將UM分為4個(gè)分子上不同、臨床上相關(guān)的亞型，這些亞型與不同的預(yù)后顯著相關(guān)[3]。盡管在尋找轉(zhuǎn)移的分子標(biāo)志物方面已經(jīng)取得了一定進(jìn)展，但UM發(fā)生遠(yuǎn)處轉(zhuǎn)移的生物學(xué)機(jī)制仍未完全明確[1-2]。

長鏈非編碼RNA(long non-coding RNA，lncRNA)是長度大于200個(gè)核苷酸的轉(zhuǎn)錄本，不能編碼蛋白質(zhì)[4]。近年來研究[5]發(fā)現(xiàn)：lncRNA能與DNA、RNA和蛋白質(zhì)等細(xì)胞大分子相互作用，參與腫瘤的發(fā)生發(fā)展。LncRNA可作為競爭性內(nèi)源性RNA(competing endogenous RNA，ceRNA)吸附microRNA(miR)，間接調(diào)控miR靶點(diǎn)mRNA的翻譯和蛋白質(zhì)合成[6]。通過這種調(diào)控方式，lncRNA促進(jìn)胃癌[7]和胰腺癌[8]等多種惡性腫瘤的發(fā)生和轉(zhuǎn)移。在人類UM組織中，已報(bào)道了幾種致癌性lncRNA(例如HOXA11-AS[9]，RHPN1-AS1[10]和PVT1[11])的顯著高表達(dá)，體外實(shí)驗(yàn)證明它們能促進(jìn)腫瘤細(xì)胞的增殖、遷移和侵襲，但尚無研究對lncRNA表達(dá)異常及其在UM轉(zhuǎn)移中的作用進(jìn)行系統(tǒng)性分析，而深入了解lncRNA在其中的作用機(jī)制可能對揭示UM轉(zhuǎn)移的潛在機(jī)制具有重要意義。

本研究的目的是利用生物信息學(xué)方法，全面解析TCGA-UVM數(shù)據(jù)集中，lncRNA的表達(dá)與相關(guān)ceRNA網(wǎng)絡(luò)在UM轉(zhuǎn)移中的作用，為進(jìn)一步實(shí)驗(yàn)研究提供生物信息學(xué)證據(jù)。

1 資料與方法

1.1 數(shù)據(jù)采集與處理

本研究從TCGA數(shù)據(jù)門戶網(wǎng)站(https://tcgadata.nci.nih.gov/tcga/)獲得所有UM患者(80例)的RNA測序數(shù)據(jù)(第3級)、miR測序數(shù)據(jù)(第3級)和臨床數(shù)據(jù)。數(shù)據(jù)使用R/Bioconductor軟件包的TCGAbiolinks程序進(jìn)行預(yù)處理，并用biomaRt包參考GENCODE(GRCh38)進(jìn)行基因注釋，以下轉(zhuǎn)錄本類型定義為lncRNA：lincRNA、antisense、sense_intronic、processed_transcript、sense_over lap ping、3prime_over lap ping_ncRNA 和macro_lncRNA。本研究最終獲得19 249 個(gè)mRNA、3 742個(gè)lncRNA和1 881個(gè)miR。根據(jù)隨訪期間是否發(fā)生遠(yuǎn)處轉(zhuǎn)移，80例UM患者分為轉(zhuǎn)移組和非轉(zhuǎn)移組兩組。

1.2 DEmRNA、DEmiR 和DElncRNA

使用R/Bioconductor包edgeR程序分析轉(zhuǎn)移組與非轉(zhuǎn)移組腫瘤樣本中DE的DEmRNA、DElncRNA和DEmiR。首先利用edgeR程序內(nèi)置的濾過功能刪除在所有樣品低表達(dá)RNA類型，再采用Benjamini-Hochberg方法計(jì)算錯誤發(fā)現(xiàn)率(false discovery rate，F(xiàn)DR)和DE倍數(shù)變化(fold change，F(xiàn)C)，將滿足FDR ＜0.05且log2(FC) ＞1.0的RNA定義為DE的RNA(即DEmRNA、DElncRNA和DEmiR)，用pheatmap工具包繪制熱圖。熱圖聚類分析采用Ward.D層次聚類法，以歐式距離定義層次距離。

1.3 ceRNA 網(wǎng)絡(luò)構(gòu)建

每一個(gè)lncRNA-miR-mRNA相互作用組合定義為一個(gè)ceRNA 單元。我們首先從miRcode 和TargetScan數(shù)據(jù)庫檢索miR與mRNA的所有調(diào)控關(guān)系，并取這兩個(gè)數(shù)據(jù)庫的交集作為最終的miRmRNA調(diào)控作用單元；從miRcode數(shù)據(jù)庫檢索miR和lncRNA的調(diào)控關(guān)系。以上數(shù)據(jù)去除重復(fù)條目后，最終獲得623622對miR-mRNA作用單元和 114 273對miR-lncRNA作用單元。根據(jù)ceRNA的作用原理，lncRNA和mRNA的表達(dá)水平應(yīng)存在正相關(guān)，因此我們計(jì)算了mRNA和lncRNA表達(dá)水平的Pearson相關(guān)系數(shù)(Pearson correlation coefficient，PCC)，并選擇PCC在最高5%百分位以內(nèi)的RNA入選構(gòu)建ceRNA網(wǎng)絡(luò)。

1.4 ceRNA 網(wǎng)絡(luò)的拓?fù)浞治龊妥泳W(wǎng)絡(luò)構(gòu)建

采用Cytoscape v3.6.0軟件構(gòu)建ceRNA網(wǎng)絡(luò)，利用軟件自帶程序分析每個(gè)節(jié)點(diǎn)(即差異RNA)的中間中心度(betweenness centrality，BC)，作為該節(jié)點(diǎn)在網(wǎng)絡(luò)中的中心度指標(biāo)，將所有節(jié)點(diǎn)的BC值排序，將BC值較高的lncRNA作為中樞性lncRNA，以這些lncRNA為核心，利用Cytoscape v3.6.0軟件再次構(gòu)建ceRNA子網(wǎng)絡(luò)。

1.5 基因功能與富集分析

使用在線工具KOBAS 3.0[12](http://kobas.cbi.pku.edu.cn)中的基因富集模塊對ceRNA網(wǎng)絡(luò)中的mRNA進(jìn)行基因功能與富集分析，將所有滿足校正P＜0.05(Corrected P Value)相關(guān)的基因本體(Gene Ontology，GO)和通路(Kyoto Encyclopedia of Genes and Genomes，KEGG)納入進(jìn)一步分析和展示，P值校正采用Benjamini-Hochberg法。

1.6 統(tǒng)計(jì)學(xué)處理

采用R軟件(3.6.0版)分析數(shù)據(jù)，連續(xù)性數(shù)據(jù)采用均數(shù)±標(biāo)準(zhǔn)差(±s)表示。將DE的RNA按其表達(dá)水平中位數(shù)分為高表達(dá)組和低表達(dá)組，應(yīng)用Kaplan-Meier生存曲線進(jìn)行生存分析。P＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 DE 的RNA 的篩選

從TCGA數(shù)據(jù)庫下載和提取到80例UM患者的數(shù)據(jù)，按隨訪期間是否發(fā)生轉(zhuǎn)移分為轉(zhuǎn)移組23例(28.8%)和非轉(zhuǎn)移組57例(71.2%)，利用edgeR算法排除低表達(dá)RNA后，篩選得到15 135個(gè)mRNA，3 051個(gè)lncRNA和719個(gè)miR進(jìn)行進(jìn)一步分析。edgeR算法篩選得到轉(zhuǎn)移組916個(gè)DEmRNA(上調(diào)346個(gè)，下調(diào)570個(gè))、72個(gè)DEmiR(上調(diào)27個(gè)，下調(diào)45個(gè))和272個(gè)DElncRNA(上調(diào)118個(gè)，下調(diào)154個(gè)；圖1)。鑒于ceRNA的作用模式，lncRNA通過間接正向調(diào)控mRNA發(fā)揮生物學(xué)作用，在轉(zhuǎn)移過程中表達(dá)升高的基因/mRNA更可能成為未來治療的靶點(diǎn)。因此，為了分析具有致癌功能的lncRNA(onco-lncRNA)及其相關(guān)ceRNA網(wǎng)絡(luò)，僅納入上調(diào)的lncRNA和mRNA行后續(xù)研究。針對DE的RNA繪制熱圖，并使用Ward.D算法和歐幾里得距離進(jìn)行聚類分析(圖1D，圖2A)。

圖1 火山圖和熱圖顯示UM轉(zhuǎn)移組差異表達(dá)RNAFigure 1 Volcano plots and heatmap showing the differentially expressed RNAs

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

為了構(gòu)建與UM轉(zhuǎn)移相關(guān)的ceRNA網(wǎng)絡(luò)，首先從miR靶點(diǎn)數(shù)據(jù)庫(即miRcode和TargetScan)獲得了405個(gè)DEmiR-DEmRNA互作單元和146個(gè)DEmiRDElncRNA互作單元，并采用PCC閾值保留mRNA和lncRNA表達(dá)顯著正相關(guān)的對子后，最終獲得67個(gè)mRNA，7個(gè)miR(miR-221和miR-222具有相同靶點(diǎn))和30個(gè)lncRNA納入ceRNA網(wǎng)絡(luò)，并形成616個(gè)lncRNA-miR-mRNA單元。所有DE的RNA組成了一個(gè)相互連通的網(wǎng)絡(luò)(圖2B)，該網(wǎng)絡(luò)包含103個(gè)節(jié)點(diǎn)和181條邊。

圖2 熱圖顯示構(gòu)成ceRNA網(wǎng)絡(luò)的差異表達(dá)RNA和ceRNA網(wǎng)絡(luò)Figure 2 Heatmap of differentially expressed RNAs in the ceRNA triples and the global view of the entire ceRNA network

2.3 GO 分析與KEGG 通路分析

使用在線工具KOBAS(v3.0)對所有上調(diào)表達(dá)的mRNA進(jìn)行GO生物過程富集分析和KEGG通路分析，共富集到374個(gè)GO類別，其中包括53個(gè)分子功能(molecular f unctions，MF)，37個(gè)細(xì)胞成分(cellular components，CC)和284個(gè)生物過程(biological process，BP)子類，其中，最相關(guān)的生物學(xué)過程包括：多細(xì)胞生物過程、解剖結(jié)構(gòu)形態(tài)發(fā)生、細(xì)胞分化調(diào)控和細(xì)胞粘附。KEGG通路分析顯示有14條通路顯著相關(guān)，包括細(xì)胞外基質(zhì)(extracellular matrix，ECM)-受體相互作用、PI3K-Akt信號通路和R ap1信號通路，圖3 顯示了P值最小的前20 個(gè)GO 類別和前10 個(gè)KEGG通路。

圖3 ceRNA網(wǎng)絡(luò)中上調(diào)mRNA的基因富集分析Figure 3 Gene enrichment analysis of all upregulated mRNAs in the ceRNA network

2.4 拓?fù)浞治龊蚦eRNA 子網(wǎng)絡(luò)分析

采用Cytoscape(v3.6)軟件中的NetworkAnalyzer模塊分析了整個(gè)ceRNA網(wǎng)絡(luò)每個(gè)節(jié)點(diǎn)的拓?fù)涮卣鳎⒂?jì)算其BC值，BC值較高的節(jié)點(diǎn)可能與UM轉(zhuǎn)移相關(guān)性更強(qiáng)，從而定義為中心RNA。在ceRNA網(wǎng)絡(luò)的103個(gè)節(jié)點(diǎn)中，BC值最高的前15個(gè)節(jié)點(diǎn)如圖4A所示，包括6個(gè)lncRNA(LINC00861、LINC02421、BHLHE40-AS1、LINC01252、LINC00513和LINC02389)和3個(gè)mRNA(UNC5D、BCL11B和MTDH)然后，通過分析與核心lncRNA直接作用的一級miR，及二級mRNA來構(gòu)建以該核心lncRNA為中心的ceRNA子網(wǎng)絡(luò)(圖4)。

圖4 ceRNA網(wǎng)絡(luò)中介中心性最高的15個(gè)節(jié)點(diǎn)及子網(wǎng)絡(luò)分析Figure 4 Top 15 nodes with the highest betweenness centrality (BC) and the subnetwork analysis of hub lncRNAs

2.5 差異RNA 表達(dá)水平與生存

Kaplan-Meier分析顯示6個(gè)核心lncRNA (LINC00861，LINC02421，BHLHE40-AS1，LINC01252，LINC00513和LINC02389)、5個(gè)miR(miR-221，miR-222，miR-506，miR-507，miR-876(除了miR-182和miR-183)和3個(gè)核心mRNA(UNC5D，BCL11B和MTDH)的表達(dá)水平與總體生存率顯著相關(guān)(圖5，6)，其中miR高表達(dá)與患者生存預(yù)后較好顯著相關(guān)(圖5)，而lncRNA和mRNA高表達(dá)與生存預(yù)后差顯著相關(guān)(圖6)。

圖5 Kaplan-Meier分析差異表達(dá)的miR與患者生存的關(guān)系Figure 5 Kaplan-Meier plot of the relationship between seven differently expressed microRNAs and the survival time of patients

圖6 Kaplan-Meier分析差異表達(dá)的核心lncRNA和核心mRNA與患者生存的關(guān)系Figure 6 Kaplan-Meier plot of the relationship between six hub lncRNAs,three hub mRNAs and the survival time of patients

3 討論

采用成熟的生物信息學(xué)算法及高質(zhì)量的TCGA-UM數(shù)據(jù)集，本研究發(fā)現(xiàn)多個(gè)與UM遠(yuǎn)處轉(zhuǎn)移相關(guān)的編碼和非編碼RNA，并以lncRNA為核心建立了一個(gè)相互聯(lián)系的ceRNA網(wǎng)絡(luò)，該網(wǎng)絡(luò)中的核心mRNA與腫瘤發(fā)生和轉(zhuǎn)移的多個(gè)GO生物學(xué)過程和KEGG通路顯著相關(guān)，核心lncRNA形成各自的ceRNA子網(wǎng)絡(luò)。

LncRNAs如何在惡性腫瘤發(fā)生和轉(zhuǎn)移中發(fā)揮調(diào)控功能的確切機(jī)制仍不清楚，體內(nèi)和體外實(shí)驗(yàn)證據(jù)逐漸證明：lncRNA可以通過吸附和抑制miR的功能，間接調(diào)節(jié)癌基因的表達(dá)，從而形成了lncRNA-miR-mRNA的競爭性內(nèi)源RNA的ceRNA新機(jī)制。迄今為止，許多研究發(fā)現(xiàn)人類惡性腫瘤中存在ceRNA調(diào)控機(jī)制[12-14]，但這種機(jī)制在UM轉(zhuǎn)移中是否存在，還缺少系統(tǒng)全面研究。本研究全面分析了TCGA-UM數(shù)據(jù)集的80個(gè)腫瘤樣本的多種類型RNA表達(dá)譜，并鑒定了570個(gè)mRNA，45個(gè)miR和154個(gè)lncRNA可能參與UM轉(zhuǎn)移，基于ceRNA作用原理和開源數(shù)據(jù)庫(miRcode和TargetScan)，最終鑒定出67個(gè)mRNA、7個(gè)miR和30個(gè)lncRNA形成了ceRNA網(wǎng)絡(luò)，這是目前第一個(gè)關(guān)于UM轉(zhuǎn)移相關(guān)ceRNA網(wǎng)絡(luò)構(gòu)建的研究。

在該ceRNA網(wǎng)絡(luò)中，miR-182、miR-183、miR-221、miR-222、miR-506、miR-507和miR-876共7個(gè)下調(diào)的miR構(gòu)成了ceRNA網(wǎng)絡(luò)的核心。根據(jù)以往研究，其中miR-506、miR-507和miR-876的表達(dá)水平與惡性腫瘤的轉(zhuǎn)移呈負(fù)相關(guān)，起抑癌基因的作用。如miR-506和miR-507可以抑制乳腺癌[15]和皮膚黑色素瘤[16]的增殖和轉(zhuǎn)移。同樣，miR-876可抑制頭頸部鱗狀細(xì)胞癌的侵襲生長[17]。這些證據(jù)與本研究在轉(zhuǎn)移性UM中發(fā)現(xiàn)的miR-506和miR-507表達(dá)降低的結(jié)果一致。但是既往研究發(fā)現(xiàn)miR-221/222在肝細(xì)胞癌[18]和宮頸鱗狀細(xì)胞癌[19]中發(fā)揮致癌作用，與本研究發(fā)現(xiàn)的轉(zhuǎn)移性UM中miR-221/222的下調(diào)相矛盾，提示miR-221/222致癌和抑癌作用的發(fā)揮可能具有一定的組織和細(xì)胞特異性。

通過ceRNA網(wǎng)絡(luò)拓?fù)涮卣鞣治?，本研究發(fā)現(xiàn)6個(gè)lncRNA(LINC00861、LINC02421、BHLHE40-AS1、LINC01252、LINC00513和LINC02389)和3個(gè)mRNA(UNC5D、BCL11B和MTDH)作為ceRNA網(wǎng)絡(luò)的中樞RNA。關(guān)于其中6種lncRNA的生物學(xué)功能的文獻(xiàn)報(bào)道較少，有研究發(fā)現(xiàn)BHLHE40-AS1[20]和LINC00861[21]的表達(dá)水平升高與乳腺癌的侵襲和肝細(xì)胞癌的早期復(fù)發(fā)有關(guān)，提示它們的致癌作用，與本研究的結(jié)果一致，但對于其他4種lncRNA(LINC02421、LINC01252、LINC00513和LINC02389)的表達(dá)水平與腫瘤發(fā)生發(fā)展的關(guān)系尚無足夠證據(jù)，它們是否與UM細(xì)胞的高度侵襲相關(guān)是將來實(shí)驗(yàn)研究的有趣方向。

與核心lncRNA 的研究相比，3 個(gè)核心mRNA(UNC5D、BCL11B和MTDH)在惡性腫瘤中的作用研究較多。有研究表明它們與多種癌癥的高轉(zhuǎn)移風(fēng)險(xiǎn)相關(guān)。如MTDH既可以促進(jìn)多發(fā)性骨髓瘤細(xì)胞增殖，又可以抑制其凋亡[22]，在腎癌[23]和結(jié)直腸癌[24]中，MTDH的高表達(dá)與遠(yuǎn)處轉(zhuǎn)移風(fēng)險(xiǎn)呈正相關(guān)。在皮膚黑色素瘤中，基因敲除MTDH可顯著抑制裸鼠黑色素瘤生長，也提示它在黑色素瘤中具有致癌作用，MTDH也可能成為治療皮膚黑色素瘤的潛在靶點(diǎn)[25]。BLC11B也是一種癌基因，在T細(xì)胞急性淋巴細(xì)胞白血病[26]和神經(jīng)膠質(zhì)瘤[27]中表達(dá)上調(diào)。這些既往報(bào)道的數(shù)據(jù)均與本研究的結(jié)果一致，但它們是否同樣參與了細(xì)胞增殖、遷移等生物學(xué)過程，還需要體內(nèi)和體外實(shí)驗(yàn)驗(yàn)證。

綜上所述，基于TCGA數(shù)據(jù)庫，本研究利用生物信息學(xué)方法鑒定出UM轉(zhuǎn)移相關(guān)的多個(gè)lncRNA，并以此構(gòu)建了ceRNA 調(diào)控網(wǎng)絡(luò)，其中6 個(gè)中樞l(wèi)ncRNA與UM總體生存率顯著相關(guān)，將來仍需要實(shí)驗(yàn)揭示它們在UM轉(zhuǎn)移中的具體作用機(jī)制。

開放獲取聲明

本文適用于知識共享許可協(xié)議(Creative Commons)，允許第三方用戶按照署名(BY)-非商業(yè)性使用(NC)-禁止演繹(ND)(CCBY-NC-ND)的方式共享，即允許第三方對本刊發(fā)表的文章進(jìn)行復(fù)制、發(fā)行、展覽、表演、放映、廣播或通過信息網(wǎng)絡(luò)向公眾傳播，但在這些過程中必須保留作者署名、僅限于非商業(yè)性目的、不得進(jìn)行演繹創(chuàng)作。詳情請?jiān)L問：https://creativecommons.org/licenses/by-nc-nd/4.0/。