甲狀腺癌非編碼RNA與mRNA差異表達(dá)譜及競(jìng)爭(zhēng)性內(nèi)源RNA調(diào)控網(wǎng)絡(luò)分析

湯喜，張成瑤，周曉紅，輦偉奇，龔靖淋，張玉蓮，周迎春，李真華

（重慶大學(xué)附屬腫瘤醫(yī)院/重慶市腫瘤研究所/重慶市腫瘤醫(yī)院 1. 頭頸腫瘤中心 2. 腫瘤轉(zhuǎn)移與個(gè)體化診治轉(zhuǎn)化研究重慶市重點(diǎn)實(shí)驗(yàn)室，重慶 400030）

甲狀腺癌（thyroid cancer）是女性第四大常見(jiàn)的惡性腫瘤[1-2]，在2017年，全球范圍內(nèi)有56870例新發(fā)甲狀腺癌病例，占所有新發(fā)腫瘤的3.4%[3]。雖然大部分甲狀腺癌傾向于生物學(xué)及臨床的惰性，擁有較好的預(yù)后，但仍有少部分甲狀腺癌預(yù)后差，生存期短以及出現(xiàn)局部及全身轉(zhuǎn)移等，嚴(yán)重威脅人類生命健康[4-6]。由于疾病過(guò)程交錯(cuò)復(fù)雜，有眾多分子參與其中，臨床中仍然難以尋找到有效的生物預(yù)測(cè)因子。侵襲和轉(zhuǎn)移是癌癥的一個(gè)標(biāo)志，主要是導(dǎo)致甲狀腺癌的復(fù)發(fā)和預(yù)后不良而導(dǎo)致死亡。對(duì)于轉(zhuǎn)移及復(fù)發(fā)性甲狀腺癌患者，目前仍缺乏有效治療手段。因此，深入了解甲狀腺癌疾病發(fā)生和發(fā)展的分子機(jī)制，探索新的潛在治療靶點(diǎn)，將有助于提高患者的生存水平。本研究通過(guò)生物信息學(xué)方法分析甲狀腺癌發(fā)生相關(guān)差異表達(dá)基因，構(gòu)建其長(zhǎng)鏈非編碼RNA（long non-coding RNA，lncRNA）、微小RNA（microRNA，miRNA）及mRNA之間的競(jìng)爭(zhēng)性內(nèi)源RNA（competing endogenous RNA，ceRNA）調(diào)控網(wǎng)絡(luò)，并結(jié)合TCGA數(shù)據(jù)庫(kù)進(jìn)行預(yù)后分析，以期獲取影響甲狀腺癌預(yù)后的關(guān)鍵非編碼RNA與mRNA及相關(guān)通路。

1 材料與方法

1.1 數(shù)據(jù)資料收集

在癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數(shù)據(jù)庫(kù)（https://cancergenome.nih.gov）中下載568例甲狀腺癌樣本并具有完整生存時(shí)間的臨床數(shù)據(jù)，使用數(shù)據(jù)轉(zhuǎn)換工具（GDC Apps提供）下載樣本的mRNAseq和miRNAseq的基因表達(dá)數(shù)據(jù)，作為后續(xù)核心差異基因的預(yù)后分析。

1.2 差異表達(dá)基因篩選

應(yīng)用R程序語(yǔ)言包（http://bioconductor.org/packages/release/bioc/html/edgeR.html）對(duì)TCGA數(shù)據(jù)庫(kù)下載的mRNAseq和miRNAseq數(shù)據(jù)信息進(jìn)行差異表達(dá)分析，篩選出甲狀腺癌組織和正常組織中差異表達(dá)的lncRNA、miRNA及mRNA。差異基因篩選標(biāo)準(zhǔn)：校正后P值＜0.01，差異倍數(shù)＞2.0。

1.3 功能富集及蛋白互作分析

采用Cytoscape v 3.5.1軟件對(duì)差異表達(dá)基因進(jìn)行GO功能富集分析參與的生物學(xué)過(guò)程及KEGG通路富集分析。采用Fisher精確檢驗(yàn)的方法，以P＜0.05及FDR＜0.05為顯著性基因富集。

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

使用Cytoscape v 3.5.1軟件構(gòu)建差異表達(dá)lncRNA、miRNA、mRNA的調(diào)控網(wǎng)絡(luò)。使用miRcode（http://www.mircode.org）預(yù)測(cè)差異表達(dá)lncRNA-miRNA關(guān)系對(duì)，使用miRTarBase（http://mirtarbase.mbc.nctu.edu.tw），miRDB（http://www.mirdb.org）和TargetScan（http://www.targetscan.org/vert_71）預(yù)測(cè)差異表達(dá)miRNA-mRNA關(guān)系對(duì)的互作關(guān)系。

1.5 差異表達(dá)基因的生存分析

將差異表達(dá)lncRNA的表達(dá)水平進(jìn)行單因素及多因素Cox回歸分析，并將差異表達(dá)lncRNA根據(jù)中位值分成高水平組和低水平組，應(yīng)用Kaplan-Meier生存曲線進(jìn)行生存分析，篩選出與甲狀腺癌總生存期顯著相關(guān)的基因，分析與基因表達(dá)與臨床參數(shù)的關(guān)系。使用受試者工作特征曲線（receiver operating characteristic curve，ROC）和ROC曲線面積（AUC），對(duì)生存分析預(yù)測(cè)的敏感性和特異性進(jìn)行評(píng)價(jià)。P＜0.05表示差異具有統(tǒng)計(jì)學(xué)意義。

2 結(jié) 果

2.1 甲狀腺癌的差異表達(dá)基因的篩選

在TCGA數(shù)據(jù)庫(kù)中下載和提取568例樣本的lncRNA數(shù)據(jù)，數(shù)據(jù)分為2個(gè)隊(duì)列：510例甲狀腺癌樣本和58例正常樣本。下載和提取573例樣本的miRNA數(shù)據(jù)，數(shù)據(jù)分為2個(gè)隊(duì)列：514例甲狀腺癌樣本和59例正常樣本。篩選出差異表達(dá)lncRNA 497個(gè)（上調(diào)93個(gè)，下調(diào)404個(gè)），差異表達(dá)miRNA 72個(gè)（上調(diào)5個(gè)，下調(diào)67個(gè)），差異表達(dá)mRNA 1097個(gè)（上調(diào)233個(gè)，下調(diào)864個(gè)）。使用R語(yǔ)言分別制作差異表達(dá)的lncRNA、miRNA、mRNA的火山圖及熱圖（圖1）。

2.2 差異表達(dá)基因的GO富集分析和KEGG通路分析

使用Cytoscape v 3.5.1軟件對(duì)1097個(gè)差異基因進(jìn)行GO生物過(guò)程富集分析和KEGG通路分析。GO生物過(guò)程富集分析結(jié)果顯示差異表達(dá)基因主要參與單組織過(guò)程、單組織細(xì)胞過(guò)程、刺激反應(yīng)等生物學(xué)過(guò)程，受體結(jié)合、分子功能調(diào)節(jié)、鈣離子調(diào)節(jié)等細(xì)胞功能，細(xì)胞、細(xì)胞膜、細(xì)胞外組件等細(xì)胞構(gòu)成過(guò)程（圖2）。KEGG通路分析結(jié)果顯示差異表達(dá)基因主要參與神經(jīng)反應(yīng)受體-配體相互作用、細(xì)胞因子-細(xì)胞因子受體相互作用、腫瘤轉(zhuǎn)錄失調(diào)等信號(hào)通路的調(diào)節(jié)（圖3）。

2.3 甲狀腺癌ceRNA調(diào)控網(wǎng)絡(luò)的構(gòu)建

使用Cytoscape v 3.5.1軟件進(jìn)行l(wèi)ncRNA-miRNA-mRNA的ceRNA互作網(wǎng)絡(luò)構(gòu)建（圖4），結(jié)果顯示有l(wèi)ncRNA-miRNA關(guān)系對(duì)72對(duì)，miRNA-mRNA關(guān)系對(duì)13對(duì)，其中共包括差異表達(dá)lncRNA 31個(gè)，差異表達(dá)miRNA 11個(gè)和差異表達(dá)mRNA 12個(gè)（圖4）。

2.4 甲狀腺癌的生存分析

在ceRNA互作網(wǎng)絡(luò)中有2個(gè)差異表達(dá)lncRNA（MIR181A2HG、OPCML-IT1），1個(gè)差異表達(dá)miRNA（miRNA-184）和2個(gè)差異表達(dá)mRNA（E2F1、SALL3）與甲狀腺癌的總體生存期明顯有關(guān)（均P＜0.05）。將甲狀腺癌樣本以中位數(shù)值為界值，分成高表達(dá)組和低表達(dá)組，采用Kaplan-Meier分析兩組患者預(yù)后相關(guān)性，結(jié)果表明低表達(dá)組具有更長(zhǎng)的生存時(shí)間（均P＜0.05）（圖5）。使用ROC曲線和AUC評(píng)價(jià)5年生存率預(yù)測(cè)的敏感性和特異性。結(jié)果顯示，差異表達(dá)的lncRNA、miRNA、mRNA的AUC值分別是0.913、0.956、0.955（圖6）。

圖1 圖1 火山圖及熱圖分析 A：差異表達(dá)lncRNA；B：差異表達(dá)miRNA；C：差異表達(dá)mRNAFigure1 XVolcano plot and heat map analysis A: Differentially expressed lncRNAs; B: Differentially expressed miRNAs; C: Differentially expressed mRNAs

圖2 GO功能富集分析Figure2 GO functional enrichment analysis

圖3 KEGG通路分析Figure3 KEGG pathway analysis

圖4 構(gòu)建LncRNA-miRNA-mRNA（ceRNA）調(diào)控網(wǎng)絡(luò)Figure4 Construction of regulatory network of lncRNA-miRNA-mRNA (ceRNA)

圖5 Kaplan-Meier生存曲線 A：差異表達(dá)lncRNA；B：差異表達(dá)miRNA；C：差異表達(dá)mRNAFigure5 Kaplan-Meier survival curves A: Differentially expressed lncRNAs; B: Differentially expressed miRNAs; C: Differentially expressed mRNAs

3 討 論

lncRNA是一類轉(zhuǎn)錄本長(zhǎng)度超過(guò)200 nt，不編碼蛋白的RNA，近年的研究[7-10]表明lncRNA位于細(xì)胞核和細(xì)胞質(zhì)中，主要通過(guò)表觀遺傳修飾、轉(zhuǎn)錄調(diào)控（轉(zhuǎn)錄后調(diào)控）和翻譯水平調(diào)控（翻譯后調(diào)控）發(fā)揮作用，也可通過(guò)與miRNA或其他細(xì)胞因子相互作用，間接影響基因表達(dá)。廣泛參與包括X染色體沉默、染色質(zhì)修飾、核內(nèi)運(yùn)輸?shù)榷喾N重要的調(diào)控過(guò)程。也有研究[11-12]表明lncRNA的表達(dá)變化與功能喪失已經(jīng)與多種惡性腫瘤在內(nèi)的人類疾病聯(lián)系在一起。miRNA是一類轉(zhuǎn)錄本長(zhǎng)度在20～25 nt的小非編碼RNA，miRNA可以抑制mRNA的翻譯或引起RNA的轉(zhuǎn)錄后翻譯的降解[13]。ceRNA假說(shuō)揭示了一種RNA間相互作用的新機(jī)制[14]。lncRNA與miRNA及其下游靶基因之間的相互調(diào)控模式與腫瘤的發(fā)生發(fā)展密切相關(guān)，成為腫瘤研究領(lǐng)域的熱點(diǎn)。miRNA作為一個(gè)轉(zhuǎn)錄后調(diào)控的重要因子，其活性可被lncRNA通過(guò)“海綿”吸附的方式調(diào)控ceRNA[15]。lncRNA作為ceRNA競(jìng)爭(zhēng)性地與miRNA結(jié)合從而影響miRNA導(dǎo)致的基因沉默，從而調(diào)節(jié)編碼基因的蛋白質(zhì)水平，參與靶基因的表達(dá)調(diào)控，在腫瘤的發(fā)生發(fā)展中發(fā)揮重要的作用。ceRNA構(gòu)建了一個(gè)復(fù)雜的作用網(wǎng)絡(luò)，正常生理狀態(tài)下ceRNA網(wǎng)絡(luò)中的各種分子之間處于一定的平衡狀態(tài)，一旦平衡被打破就會(huì)導(dǎo)致疾病的發(fā)生。目前有研究揭示了ceRNA調(diào)控網(wǎng)絡(luò)在部分腫瘤中的主要作用[16-18]。目前仍缺乏在全基因組范圍內(nèi)，特別是基于高通量檢測(cè)與大規(guī)模樣本量的甲狀腺癌相關(guān)的lncRNA和miRNA及的ceRNA的綜合分析。

MIR181A2HG（MIR181A2 host gene，MIR181A2宿主基因）是位于染色體9q33.3的lncRNA，miR-181a 從兩條人類染色體的lncRNA轉(zhuǎn)錄而來(lái)，分別是1號(hào)染色體的MIR181A1HG和9號(hào)染色體的MIR181A2HG。有報(bào)道MIR181A2HG可以促進(jìn)子宮相關(guān)炎癥通路的發(fā)生[19]。OPCMLIT1（OPCML intronic transcript 1，OPCML內(nèi)轉(zhuǎn)錄1）是位于11號(hào)染色體的lncRNA。目前尚無(wú)MIR181A2HG和OPCML-IT1在甲狀腺癌中的表達(dá)及作用的相關(guān)報(bào)道，有待于我們進(jìn)一步進(jìn)行基礎(chǔ)實(shí)驗(yàn)和臨床試驗(yàn)。miR-184近年來(lái)作為一個(gè)新的miRNA在許多腫瘤中表現(xiàn)出異常的表達(dá)，包括肝細(xì)胞癌（作為肝細(xì)胞癌的致癌調(diào)節(jié)劑）和腎細(xì)胞癌（作為腫瘤抑制因子）等[20-21]。E2F1是炎癥控制的關(guān)鍵參與者，它與成視網(wǎng)膜細(xì)胞瘤蛋白（pRb）息息相關(guān)。pRb通過(guò)對(duì)E2F家族轉(zhuǎn)錄因子的誘導(dǎo)來(lái)控制細(xì)胞周期。有研究[22]表明變異的pRb通過(guò)激活E2F1來(lái)抑制甲狀腺腫瘤。Sall3（spalt like transcription factor 3）編碼C2H2鋅指蛋白，這種蛋白質(zhì)存在于進(jìn)化保守的基因家族中，包括果蠅、新桿狀線蟲(chóng)和脊椎動(dòng)物[23]。近期有研究[24]表明SALL3的表達(dá)與致癌作用之間有關(guān)系。研究表明在肝細(xì)胞癌中SALL3被DNA甲基化所抑制，其蛋白質(zhì)與DNA甲基轉(zhuǎn)移酶3α（DNMT3A）相互作用。另一個(gè)研究[25]則指出SALL3和HPV感染的異常超甲基化對(duì)宮頸癌的致癌作用有促進(jìn)作用。SALL3在甲狀腺癌中是否有促癌作用尚未見(jiàn)報(bào)道。

通過(guò)對(duì)lncRNA和miRNA的靶基因預(yù)測(cè)和功能分析發(fā)現(xiàn)，GO生物過(guò)程富集分析結(jié)果顯示差異表達(dá)基因主要參與單組織過(guò)程、單組織細(xì)胞過(guò)程、刺激反應(yīng)等生物學(xué)過(guò)程，受體結(jié)合、分子功能調(diào)節(jié)、鈣離子調(diào)節(jié)等細(xì)胞功能，細(xì)胞、細(xì)胞膜、細(xì)胞外組件等細(xì)胞構(gòu)成過(guò)程。KEGG通路分析結(jié)果顯示差異表達(dá)基因主要參與神經(jīng)反應(yīng)受體-配體相互作用、細(xì)胞因子-細(xì)胞因子受體相互作用、腫瘤轉(zhuǎn)錄失調(diào)等信號(hào)通路的調(diào)節(jié)。因此，初步探討與預(yù)后相關(guān)的lncRNA和miRNA的靶基因生物過(guò)程富集和信號(hào)通路分析，對(duì)于揭示甲狀腺癌發(fā)生和發(fā)展的分子機(jī)制以及腫瘤的診斷治療具有一定的臨床意義。

本研究對(duì)T C G A數(shù)據(jù)庫(kù)中甲狀腺癌的lncRNA，miRNA和mRNA測(cè)序數(shù)據(jù)與甲狀腺癌的臨床信息進(jìn)行了相關(guān)的統(tǒng)計(jì)學(xué)分析。其優(yōu)勢(shì)在于研究樣本量大，測(cè)試平臺(tái)統(tǒng)一，結(jié)果全面可靠。但其結(jié)果具有一定局限性，需要在基礎(chǔ)實(shí)驗(yàn)和臨床試驗(yàn)中進(jìn)一步驗(yàn)證。

綜上所述，通過(guò)對(duì)甲狀腺癌患者的lncRNA，miRNA和mRNA數(shù)據(jù)進(jìn)行分析發(fā)現(xiàn)，差異表達(dá)lncRNA（MIR181A2HG、OPCML-IT1），差異表達(dá)miRNA（miRNA-184）和差異表達(dá)mRNA（E2F1、SALL3）可能在甲狀腺癌的發(fā)生發(fā)展中發(fā)揮作用，具有成為新的預(yù)測(cè)甲狀腺癌預(yù)后分子標(biāo)志物的潛能。對(duì)于揭示甲狀腺癌發(fā)生、發(fā)展的分子機(jī)制及腫瘤的診斷和分子靶向治療具有一定的臨床意義。