利用TCGA數據庫構建腎透明細胞癌相關miRNA預后模型

高 艾，王昕苑，蘇依琳，蘇龍龍，張建輝，牛曉辰

(山西醫科大學,太原 030000)

腎細胞癌(簡稱腎癌)是一類常見的惡性腫瘤，起源于腎小管上皮，約占成人腫瘤發病率的2%～3%[1]，主要病理類型有腎透明細胞癌、腎乳頭狀癌、腎嫌色細胞癌以及集合管癌等[2]。其中，腎透明細胞癌(Kidney renal clear cell carcinoma,KIRC)占腎癌的80%～90%，是最主要的病理類型[3]。隨著腫瘤檢測手段與治療方法的進步，腎透明細胞癌患者的存活率得到顯著提高，但總體生存期以及無進展生存期仍然較低，超過1/3的病人在診斷時已經發生轉移，因此迫切需要新的早期診斷標志物以及治療靶點來改善患者的生存預后情況[4]。微小RNA(microRNA，miRNA)是一類廣泛存在于動植物體內的短鏈非編碼RNA，僅有21～23個核苷酸序列，通過抑制翻譯或促進細胞質中mRNA降解來調控轉錄后的基因表達[5]。大量研究證實miRNA在癌癥的發生、發展、轉移和患者生存中發揮重要作用，作為一種可能的非侵入性生物標志物，其預后價值在多種腫瘤中被廣泛研究[6-8]。但相關研究報道在腎透明細胞癌中較少，僅有的幾個研究也只是從miRNA的角度進行預后評估[9-10]，未能聯系下游mRNA，進而繪制miRNA-mRNA調控網絡，以miRNA-mRNA關系對的形式深入挖掘其在腎透明細胞癌中的作用，為相關診斷與治療提供參考，而這恰恰成為本研究的主要目的。

1 材料與方法

1.1 數據獲取與整理

于2020年7月12日檢索癌癥基因組圖譜數據庫(The Cancer Genome Atlas,TCGA,https://portal.gdc.cancer.gov/)并下載KIRC的mRNA與miRNA測序信息和相應臨床數據，其中包括539例KIRC患者的腫瘤組織測序信息和72例正常組織測序信息。利用perl語言(perl 5.30.2,http://www.perl.org/)分別將mRNA與miRNA測序數據合并為矩陣文件；根據Ensembl數據庫(http://www.ensembl.org/)將mRNA文件中“Ensembl_Stable_ID”轉換為“Gene Symbol”，后整理為獨立文件；根據miRbase數據庫(http://www.mirbase.org/)將miRNA文件中的“MIMAT_ID”轉換為“hsa-miR_ID”，后整理為獨立文件。由于TCGA數據庫屬于公共的開放獲取數據，故本研究不需要相關倫理學審核與批準。

1.2 差異分析

利用R語言(R 3.6.3,https://www.r-project.org/)中的edgeR包對mRNA與miRNA分別進行差異分析，設置FDR(BH)矯正后的閾值P.adj<0.05，對數差異表達倍數變化絕對值|log2FC|>2；利用pheatmap包分別繪制表達上調或下調的前20個基因的熱圖。

1.3 預后模型構建

利用perl語言將差異表達的miRNA同生存時間與生存狀態進行合并，刪除臨床數據不完整以及生存時間小于30天的樣本；首先設置生存分析過濾條件：P.adj<0.05，后利用R語言中的caret包將樣本隨機分為實驗組(Train組)和驗證組(Test組)；利用survival包對Train組進行單因素與多因素Cox回歸分析，篩選與患者預后密切相關的miRNA并構建預后風險評分模型，模型計算公式為：風險值(risk score)=風險基因表達量1×coef1+風險基因表達量2×coef2+...+風險基因表達量n×coefn(coef為風險系數)；利用預后模型分別計算Train組與Test組各樣本的風險值。將各組樣本的風險值從低到高排序，依據中位數將患者分為低風險組和高風險組，風險值越大，生存率越低。利用survminer與survival ROC包對Train組與Test組進行Kaplan-Meier (K-M)生存分析并計算ROC曲線的AUC值(設置Train組：P.adj<0.01，AUC>0.70；Test組：P.adj<0.01，AUC>0.69)，如分組不滿足設置條件，則重新利用caret包進行分組并循環。

1.4 預后模型評價

利用survival包輸出納入模型的miRNA的生存曲線，驗證其生存情況與coef的關系；同時利用OncoLnc網頁工具(http://www.oncolnc.org/)輸出模型中miRNA的生存曲線，驗證數據準確性(截斷值：50%)；根據預后模型計算所有樣本、Train組以及Test組的風險值，由低到高排序后根據中位數將患者分為高低風險兩組，輸出各組中模型的生存曲線，驗證模型的準確性；利用survival ROC包輸出所有樣本、Train組以及Test組的ROC曲線，評價模型的預測能力。將所有樣本的生存時間、生存狀態、年齡、性別、WHO分級、分期、TNM分期及風險評分整理為數值數據，其中：生存狀態(0，存活；1，死亡)、性別(0，女性；1，男性)，利用survival包進行單因素與多因素的獨立預后分析，從而評估風險值是否能作為患者的獨立預后因子(評價標準：P.adj<0.01)。

1.5 miRNA-mRNA調控網絡構建

利用miRDB(http://mirdb.org/)、miRTarBase(http://mirtarbase.cuhk.edu.cn/)與TargetScan(http://www.targetscan.org/)數據庫對納入模型的miRNA進行靶基因預測，篩選同時被兩個及以上數據庫收錄的靶基因與1.2中差異表達的mRNA取交集，明確其靶關系以及上下調表達關系，利用Cytoscape 3.7.2軟件(https://cytoscape.org/)繪制miRNA-mRNA調控網絡。利用perl語言將納入網絡的mRNA同臨床數據進行合并，利用R語言中的survival包進行生存分析(設置P.adj<0.05)，篩選與患者生存預后密切相關的mRNA。為了進一步篩選與患者生存預后相關的miRNA-mRNA調控關系對，根據miRNA對mRNA的功能發揮抑制作用這一生物學基礎，制定關系對篩選標準：“高表達，生存率低”的miRNA對應的mRNA應符合“低表達，生存率低”；“低表達，生存率低”的miRNA對應的mRNA應符合“高表達，生存率低”。這樣的miRNA-mRNA調控關系對具有重要生物學意義，通過選擇性調控上游miRNA的表達，進而影響mRNA的表達，可改善腎透明細胞癌患者的生存情況。

2 結果

2.1 差異分析

相比于正常組織，在KIRC腫瘤組織中共有3 613個差異表達的mRNA(上調2 603個，下調1 010個)，上調與下調最顯著的前20個mRNA(見圖1a)(C代表正常對照組，T代表腫瘤組)；共有49個差異表達的miRNA(上調28個，下調21個)，上調與下調最顯著的前20個miRNA(見圖1b)。

圖1 差異表達的mRNA與miRNA(C代表正常對照組，T代表腫瘤組)Fig.1 Differentially expressed mRNAs and miRNAs(C represents control group,T represents tumor group)

2.2 預后模型構建

刪除臨床數據不完整以及生存時間小于30天的樣本，共有499個樣本的miRNA表達數據被納入分析。單因素Cox回歸分析顯示，共有3個miRNA的表達情況與患者的生存預后明顯相關(見表1)；多因素Cox回歸分析顯示，共有2個miRNA被納入預后風險評分模型(見表2)，分別是hsa-miR-21-5p和hsa-miR-1251-5p；預后風險評分模型計算公式為：風險值(risk score)=hsa-miR-21-5p表達量×0.603+hsa-miR-1251-5p表達量×-0.093。模型構建所基于的Train組共有251個樣本，Test組共有248個樣本。

表1 單因素Cox回歸分析Table 1 Univariate Cox regression analysis

表2 多因素Cox回歸分析Table 2 Multivariate Cox regression analysis

2.3 預后模型評價

hsa-miR-21-5p生存曲線(見圖2a)顯示，高表達組，生存率低；hsa-miR-1251-5p生存曲線(見圖2b)顯示，低表達組，生存率低。根據預后模型計算所有樣本、Train組以及Test組的風險值，所有樣本(見圖2c)、Train組(見圖2d)以及Test組(見圖2e)的生存曲線均顯示，高風險組生存率明顯低于低風險組，模型準確。所有樣本(見圖2f)、Train組(見圖2g)以及Test組(見圖2h)的ROC曲線顯示，模型具備一定預測能力。利用OncoLnc網頁工具輸出hsa-miR-21-5p(見圖3a)與hsa-miR-1251-5p(見圖3b)在KIRC中的生存曲線，結果與本研究一致。單因素(見圖4a)與多因素(見圖4b)獨立預后分析顯示，年齡與風險值均可作為患者生存預后的獨立預測因子，且隨著二者的增大，患者生存率逐漸降低。

圖2 生存曲線與ROC曲線Fig.2 Survival curves and ROC curves

圖3 OncoLnc網頁工具Fig.3 OncoLnc web tool

圖4 單因素與多因素獨立預后分析Fig.4 Univariate and multivariate independent prognostic analyses

2.4 miRNA-mRNA調控網絡構建

利用miRDB、miRTarBase與TargetScan三個數據庫對納入模型的2個miRNA，hsa-miR-21-5p與hsa-miR-1251-5p進行靶基因預測，共預測到12 354個靶基因，其中同時被兩個及以上數據庫收錄的靶基因有991個，與1.2中差異表達的mRNA取交集，共得到31個靶基因，所構建miRNA-mRNA調控網絡(見圖5)(圖中紅色代表表達上調，綠色代表表達下調)。對網絡中的mRNA進行生存分析，共有8個mRNA的表達與患者生存預后相關，分別是：AIF1L、ALX1、CLIC5、COBLL1、FCGR1A、FCGR1B、FGF1、NIPAL1，生存曲線(見圖6)。根據制定的miRNA-mRNA關系對篩選標準，共有7個關系對具有重要生物學意義，通過調控miRNA的表達，可影響下游mRNA水平的變化，進而改善患者的生存預后情況(見表3)。

圖5 miRNA-mRNA調控網絡Fig.5 miRNA-mRNA regulatory network

圖6 mRNA生存曲線Fig.6 Survival curves of mRNAs

表3 miRNA-mRNA調控關系對Table 3 miRNA-mRNA regulatory pairs

3 討 論

miRNA作為一類具有調控功能的非編碼RNA，被證明幾乎參與腫瘤發生進展的全過程，包括持續增殖、無限復制、回避生長抑制、抵抗凋亡、血管生成以及腫瘤微環境構建等，在不同類型腫瘤中表達異常的miRNA通過作用于特定的靶基因mRNA，發揮抑癌或促癌作用[11]。由于其分子量較小，結構穩定不易降解，包裹在外泌體中的miRNA極易進入循環系統，可穩定存在于血清、血漿、腦脊液、唾液以及尿液等，這些特點使其逐步成為一種可靠的生物標志物，被廣泛應用在各類腫瘤的診斷與預后評估中[12]。而最新研究表明，一種新型的細胞內活檢技術，僅需約10分鐘，即可從活細胞中分離出目標miRNA，大大縮短檢測所需時間，同時操作簡單，靈敏度和準確性都得以提升，這項技術為miRNA在體內的檢測提供便利，同樣有助于癌癥的早期診斷與評估[13]。

本研究所構建的預后模型共納入2個miRNA，分別是hsa-miR-21-5p與hsa-miR-1251-5p。miR-21是一個廣泛表達的miRNA，在多種腫瘤細胞中增高，可通過抑制PTEN或其他抑癌基因來發揮其致癌作用。研究發現，在腎透明細胞癌組織中，miR-21-5p的過表達可抑制SATB1的表達，進而導致患者的不良預后[14]；另一項研究則對石蠟包埋的腎透明細胞癌組織中miRNA的表達與患者生存情況之間的關系進行回顧性隊列研究，結果發現miR-21-5p 與 miR-210-3p表達上調最顯著，且與患者的不良預后密切相關[15]，這些研究均提示miR-21-5p可成為該腫瘤的有效預后評價指標，構建預后模型具有一定臨床意義。一項有關卵巢癌的研究發現，miR-1251-5p在癌組織中高表達，可通過靶向抑癌基因TBCC，進而促進癌細胞的增殖和自噬的發生[16]。而在本研究發現，miR-1251-5p在KIRC組織中低表達，但其高表達組生存率更高，且模型中coef值為負，即表達越多，患者的風險值越低，生存預后情況越好，提示過表達miR-1251-5p在KIRC中可能發揮抑制癌組織的作用，這與卵巢癌中的研究不一致。由于同一miRNA在不同組織中往往結合不同靶基因發揮不同功能，因此miR-1251-5p與KIRC生存預后的關系仍有待考證，具有一定研究價值。

由于一個miRNA的靶基因往往有多個，一個靶基因又可能和多個miRNA相互作用，組成復雜的調控網絡，而單純糾正某個miRNA異常可能會引起許多副作用，只有從miRNA-mRNA調控關系對的角度出發，精準調控mRNA的表達，才能實現對腫瘤的靶向治療。因此本研究通過對KIRC腫瘤組織與正常組織mRNA的表達情況進行差異分析，篩選到3 613個差異表達的mRNA，后與miRNA的靶基因進行取交集，得到31個mRNA，對這些mRNA進行生存分析并根據關系對篩選標準，篩選到7個miRNA-mRNA進行后續研究。

其中，與miR-21-5p結合的mRNA共有5個，分別是AIF1L、CLIC5、COBLL1、FGF1、NIPAL1，且在KIRC組織中均低表達，生存分析均顯示低表達組的生存率較低，提示其均可能發揮抑癌作用。(1)AIF1L(allograft inflammatory factor 1 like，同種異體移植炎癥因子1樣)是一類具有EF手型模序結構的鈣結合蛋白，其在腎臟組織中高表達。一項有關乳腺癌的研究表明，腫瘤組織中AIF1L低表達且高度甲基化，其可導致患者的不良預后，過表達AIF1L則可抑制細胞的擴散，改變細胞形態，減少突起的形成[17]。(2)CLIC5(chloride intracellular channel 5，細胞內氯離子通道5)所編碼的蛋白質與肌動蛋白為基礎的細胞骨架結構相關，同時參與腎小球足細胞和內皮細胞的維持。研究表明，在肝癌細胞Huh7中抑制CLIC5表達，可導致腫瘤細胞遷移和侵襲能力下降[18]；而在小兒急性淋巴細胞白血病中發現，過表達CLIC5可促進氧化應激誘導的DNA損傷積累，從而促進白血病的發生[19]。(3)COBLL1(cordon-bleu WH2 repeat protein like 1，藍帶WH2重復蛋白1)在腎臟組織中高表達，其具有肌動蛋白結合結構域，在臨床前列腺癌組織中表達上調，并與前列腺癌患者不良預后有關[20]；在慢性髓系白血病的急性轉化期，COBLL1高表達，且與生存率降低顯著相關[21]。(4)FGF1(fibroblast growth factor 1，成纖維細胞生長因子1)是成纖維細胞生長因子家族的一員，具有廣泛的有絲分裂和細胞存活活性，參與多種生物過程，包括胚胎發育、細胞生長、形態形成、組織修復、腫瘤生長和侵襲等。研究證實，FGF家族成員在功能上參與了腎透明細胞癌的進展，通過自分泌和(或)旁分泌的作用直接刺激細胞增殖[22]。(5)NIPAL1(NIPA like domain containing 1，NIPA類域包含子1)在食管癌組織中表達降低，過表達該基因通過阻滯細胞G2/M期，抑制細胞有絲分裂及減少細胞骨架中絲狀偽足和板狀偽足的形成，抑制癌細胞的增殖和遷移[23]；而在口腔鱗狀細胞癌中的研究則發現，NIPAL1可調節腫瘤細胞和內皮細胞的生長和粘附，發揮促癌作用[24]。因此，AIF1L極有可能是KIRC的一個抑癌基因，成為有效治療靶點；NIPAL1在不同腫瘤類型中發揮促癌或抑癌作用，其與KIRC的關系有待進一步研究與考證；而CLIC5、COBLL1、FGF1三個mRNA在幾類癌癥中均發揮促癌作用，與我們的假設相異。

與miR-1251-5p結合的mRNA有2個，分別是FCGR1A(Fc fragment of IgG receptor Ia，免疫球蛋白G Fc段受體1α)與FCGR1B(Fc fragment of IgG receptor Iβ，免疫球蛋白G Fc段受體1β)，它們是位于1號染色體上相互關聯的兩個基因，該家族共有三個亞型，共同編碼高親和力的Fc -γ受體。二者在KIRC組織中均高表達，且生存分析顯示高表達組生存率較低，提示其可能發揮促癌作用。一項前瞻性研究發現，FCGR1A的高表達與三陰性乳腺癌患者的不良預后有關[25]；另一項研究則發現，抗癌基因miR-29可通過靶向FCGR1B，降低KIRC腫瘤細胞的遷移和侵襲能力，改善患者的生存預后情況[26]。

本研究發現一個新的miRNA與KIRC生存預后相關，即hsa-miR-1251-5p；發現7個miRNA-mRNA調控關系對可為腎透明細胞癌的研究與治療提供靶點，其中miR-21-5p-AIF1L、miR-21-5p-NIPAL1、miR-1251-5p-FCGR1A與miR-1251-5p-FCGR1B，此4個關系對均未在KIRC中研究報道，具有較大研究價值與可信度。本研究的特色之處在于，利用miRNA進行預后模型構建，作為一種無創性檢查手段具有一定應用價值；所分析的miRNA與mRNA數據均來源于同一批患者，具有一定說服力；在miRNA-mRNA調控網絡的基礎上，進一步進行mRNA的生存分析，篩選到更具研究價值的調控關系對。本研究的不足之處在于，未在其他數據集中對預后模型進行驗證與評估；所篩選到的mRNA表達與生存情況未在蛋白質層面進行分析；未對有研究價值的調控關系對進行功能實驗驗證；這些將在后續研究中進一步完善。

4 結 論

1)通過對TCGA數據庫中KIRC的mRNA與miRNA數據進行差異分析，結合臨床數據通過miRNA表達量有關的單因素與多因素Cox回歸分析，成功構建可預測患者生存預后情況的風險評分模型，計算公式為：風險值(Risk score)=hsa-miR-21-5p表達量×0.603+hsa-miR-1251-5p表達量×-0.093。

2)通過對納入模型的miRNA進行靶基因預測，結果與差異表達的mRNA取交集，成功構建miRNA-mRNA調控網絡，對網絡中的mRNA進行生存分析，篩選到7個重要的調控關系對，可為相關研究與治療提供理論支持。