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

（重慶大學附屬腫瘤醫(yī)院/重慶市腫瘤研究所/重慶市腫瘤醫(yī)院 1. 頭頸腫瘤中心 2. 腫瘤轉移與個體化診治轉化研究重慶市重點實驗室，重慶 400030）

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

1 材料與方法

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

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

1.2 差異表達基因篩選

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

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

采用Cytoscape v 3.5.1軟件對差異表達基因進行GO功能富集分析參與的生物學過程及KEGG通路富集分析。采用Fisher精確檢驗的方法，以P＜0.05及FDR＜0.05為顯著性基因富集。

1.4 ceRNA調控網(wǎng)絡的構建

使用Cytoscape v 3.5.1軟件構建差異表達lncRNA、miRNA、mRNA的調控網(wǎng)絡。使用miRcode（http://www.mircode.org）預測差異表達lncRNA-miRNA關系對，使用miRTarBase（http://mirtarbase.mbc.nctu.edu.tw），miRDB（http://www.mirdb.org）和TargetScan（http://www.targetscan.org/vert_71）預測差異表達miRNA-mRNA關系對的互作關系。

1.5 差異表達基因的生存分析

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

2 結 果

2.1 甲狀腺癌的差異表達基因的篩選

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

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

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

2.3 甲狀腺癌ceRNA調控網(wǎng)絡的構建

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

2.4 甲狀腺癌的生存分析

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

圖1 圖1 火山圖及熱圖分析 A：差異表達lncRNA；B：差異表達miRNA；C：差異表達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 構建LncRNA-miRNA-mRNA（ceRNA）調控網(wǎng)絡Figure4 Construction of regulatory network of lncRNA-miRNA-mRNA (ceRNA)

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

3 討 論

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

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

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

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

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