m 6ARNA甲基化調(diào)節(jié)基因預測甲狀腺癌的臨床預后的作用

吳洪淵，王媛媛，傅 友，戴 萌，吳依芬，王 剛

南方醫(yī)科大學南方醫(yī)院1放療科，2健康管理科，廣東 廣州510515；南方醫(yī)科大學附屬東莞市人民醫(yī)院3腫瘤內(nèi)科，4放射科，廣東 東莞523059

甲狀腺癌是一種最常見的內(nèi)分泌相關(guān)的惡性腫瘤，超過90%內(nèi)分泌系統(tǒng)腫瘤都是甲狀腺癌，其發(fā)病率在逐年上升[1-5]。甲狀腺癌雖然不是一種高發(fā)腫瘤，但仍然是女性第六大腫瘤，僅次于乳腺癌、肺癌、直腸癌、宮頸癌和皮膚惡性腫瘤，女性的發(fā)病率是男性的2～3倍[6-7]。甲狀腺癌死亡率較低，其10年生存率90%～95%[8]。雖然手術(shù)、放療、左甲狀腺素治療方法在不斷進步[9-10]，仍有超過10%甲狀腺癌患者會復發(fā)[11-14]，其余治療手段對于局部進展和遠處轉(zhuǎn)移的甲狀腺癌效果有限，因此進一步明確甲狀腺癌的分子機制有助于腫瘤預后的預測及腫瘤的精準治療。

目前THCA的發(fā)病機制未完全闡述清楚。在基因水平對甲狀腺癌進行研究發(fā)現(xiàn)，RET、BRAF等基因在甲狀腺癌的發(fā)生、發(fā)展以及轉(zhuǎn)移中發(fā)揮重要作用[15]。m6A作為mRNA中最常見的甲基化修飾，無論是在mRNA還是lncRNA，都存在大量的m6A修飾。m6A對癌癥的作用反應(yīng)在對癌癥相關(guān)基因表達的調(diào)節(jié)中，可促進腫瘤細胞的生長及侵襲，亦可產(chǎn)生抑制作用。m6ARNA調(diào)節(jié)基因METTL3與乳腺癌、卵巢癌和肝癌的發(fā)生發(fā)展關(guān)系緊密[16-18]。此外，F(xiàn)TO和ALKBH5在乳腺癌中同樣發(fā)揮重要作用[19-20]。雖然越來越多的研究證明m6ARNA甲基化在眾多腫瘤的發(fā)生發(fā)展中的作用，但是其與THCA的關(guān)系仍不清楚，也沒有可依靠m6A構(gòu)建的預后模型。本研究使用TCGA數(shù)據(jù)分析了THCA與正常組織中，20個m6ARNA甲基化調(diào)節(jié)基因的表達差異，用Lasso回歸分析選出4個基因構(gòu)建風險模型。然后驗證其在THCA中的預測作用。本研究發(fā)現(xiàn)m6ARNA調(diào)節(jié)基因在THCA的發(fā)生發(fā)展中起著重要的作用，并且可作為獨立預測基因來預測THCA病人的預后，為進一步個體化治療提供預測條件。

1 資料與方法

1.1 數(shù)據(jù)下載

從TCGA數(shù)據(jù)庫下載RNA-seq表達數(shù)據(jù)（https://portal.gdc.cancer.gov/），納入58個正常組織，497個腫瘤組織?；颊吲R床信息（表1）。

1.2 m6A甲基化調(diào)節(jié)基因選擇

選擇20個被公認的m6A調(diào)節(jié)基因，從TCGA數(shù)據(jù)庫中得到表達數(shù)據(jù)，在正常組織和腫瘤組織中用limma包用wilcox檢驗方法進行差異基因篩選，然后將臨床病理資料進行比較。

1.3 生物信息學分析

所有生物信息分析都在R語言軟件中進行。熱圖繪制采用pheatmap包，用corrplot包用Pearson相關(guān)性分析對20個m6A調(diào)節(jié)基因進行相關(guān)性分析。將20個m6A調(diào)節(jié)基因用glmnet包進行l(wèi)asso風險回歸分析，選取4個與預后相關(guān)基因，并根據(jù)風險值將病人分為高低風險組，用survival包進行高低風險組生存分析及survivalROC包進行ROC分析檢驗模型的有效性，Cox回歸分析用來分析風險值與年齡、性別、stage分期、T分期、N分期之間的關(guān)系。森林圖用forestplot包進行繪制。為了進行蛋白水平驗證，進入The Human Protein Atlas（https://www.proteinatlas.org/）進 行 IGF2BP1、IGF2BP2、FTO與RBM15在正常組織和腫瘤組織中的免疫組化表達情況探索。

1.4 預后模型驗證

隨機選取原數(shù)據(jù)集中的一半患者，根據(jù)風險值將其分為高風險組和低風險組，在這兩組中進行OS及ROC分析，并進行Cox分析評估風險值與年齡、性別、stage分期、T分期、N分期之間的關(guān)系。

1.5 統(tǒng)計學分析

用wilcoxon test統(tǒng)計方法分析腫瘤樣品和正常樣品中m6ARNA甲基化調(diào)節(jié)基因的差異表達。區(qū)分病人為高低風險的cutoff值為上面得到的風險值。采用Kaplan-Meier方法分析OS，利用卡方檢驗來分析風險值和臨床特征變量。最后構(gòu)建出的風險模型在臨床亞組中進行生存分析驗證。

2 結(jié)果

2.1 THCA中m6A調(diào)節(jié)基因的表達

對20個m6A調(diào)節(jié)基因進行分析，結(jié)果顯示，17個m6A RNA調(diào)節(jié)基因都有表達差異。IGF2BP2、IGF2BP3、HNRNPC表達上調(diào)，而YTHDC2、YTHDC1、ALKBH5、HNRNPA2B1、ZC3H13、FTO、METTL14、WTAP、YTHDF1、RBM15、KIAA1429、METTL3、YTHDF3、IGF2BP1表達下調(diào)（圖1）。

2.2 m6A甲基化調(diào)節(jié)基因相關(guān)性

20個調(diào)節(jié)基因中除了IGF2BP1和IGF2BP3，剩下18個基因都具有相關(guān)性，其中相關(guān)性最大的是YTHDF3和KIAA1429（圖2）。

2.3 m6A調(diào)節(jié)基因在THCA中的預后作用

20個基因進行l(wèi)asso回歸分析后篩選得到4個基因，分別是IGF2BP2（-9.95183573118394e-05）；

FTO（0.000109707439996991）；

IGF2BP1（0.0033484012761759）；

RBM15（0.00063213881736861）（圖3A、B）。根據(jù)風險值的中位值，將病人分為高危組和低危組，高風險組OS與低風險組相比明顯較低（圖3C，P=4.278e-02）。ROC曲線用來評估預測模型的敏感性和特異性，結(jié)果顯示AUC值為0.783（圖3D）。

為驗證所篩選出的4個基因在正常組織和腫瘤組織中的蛋白水平的表達，從The Human ProteinAtlas官網(wǎng)進行四個基因檢索，結(jié)果顯示：IGF2BP1在甲狀腺正常組織和腫瘤組織中均為弱表達；IGF2BP2在甲狀腺正常組織中表達高于腫瘤組織；FTO在正常組織和腫瘤組織中均高表達；RBM15未找到相關(guān)表達情況（圖3E～J）。

2.4 THCA病人中風險值的作用

高低風險組病人臨床信息（表1）。熱圖顯示選擇出來的4個基因在高低風險組中的表達量及與臨床特征的相關(guān)性（圖4A）。結(jié)果顯示，在stage T分期（P＜0.05）和N分期（P＜0.01）中差異有統(tǒng)計學意義。單因素分析顯示，年齡（以60歲劃分為兩組）、stage、T、riskScore與OS相關(guān)（P＜0.001）；多因素分析顯示，年齡和riskScore與OS相關(guān)（P＜0.001，圖4B、C）。

為進一步研究riskScore在臨床亞組中的預測價值，根據(jù)病人臨床特征進行分組（stage、T、N、M、age），高風險和低風險組的OS在N0（P=4.911e-02），stage3/4（P=4.322e-02）亞組中有差異，且高風險組OS更低。其他分組未見明顯異常。

2.5 預測模型的驗證

隨機選中251例患者納入為驗證組，根據(jù)風險值將患者分為高風險組（n=121）、低風險組（n=130）。生存分析顯示，與低風險組相比，高風險組的OS更短（P=8.756e-03，圖6A），ROC顯示在驗證組中AUC值為0.778（圖6B）。在驗證組中，單因素分析表明年齡、stage分期、T分期和風險值與OS明顯相關(guān)（圖6C）。

圖1 m6ARNA甲基化調(diào)節(jié)基因在THCA中的表達Fig.1 The expression of m6A RNA methylation regulators in THCA The heat map(A)and The violin plot(B)of 20 m6A RNA methylation regulators between THCA tumor samples and normalsamples.*P＜0.05,**P＜0.01,***P＜0.001.

圖2 THCA中20個m6ARNA甲基化調(diào)節(jié)基因的相關(guān)性Fig.2 The correlation between 20 m6A RNA methylation regulators in THCA.

3 討論

多組學研究表明THCA是一種分子水平高度異質(zhì)性的疾病[21-22]，分子表達可作為一種THCA診斷、預后及治療反應(yīng)的潛在預測因素，可克服臨床和組織病理學特征的局限性。mRNA、miRNA和lncRNA已被廣泛探索，很多mRNA、miRNA和lncRNA相關(guān)的標志已已有研究證明可提高THCA診斷和預后[23-28]。多種基因總是相互作用來調(diào)節(jié)腫瘤的發(fā)生發(fā)展[29]。最近幾年，m6ARNA甲基化正成為不同腫瘤研究領(lǐng)域中的熱點，其與mRNA、miRNA和lncRNA之間的作用在很多腫瘤中也被廣泛研究。m6ARNA甲基化調(diào)節(jié)基因在腫瘤的發(fā)生發(fā)展中起著重要的作用[30]，但其在THCA中的作用機制不明，因此很有必要探討m6ARNA調(diào)節(jié)基因?qū)HCA的作用。

圖3 THCA中20個m6ARNA甲基化調(diào)節(jié)基因的風險特征Fig.3 The risk signature of 20 m6ARNA methylation regulators in THCA.A,B:The coefficients calculated by Lasso Cox regression analysis;C:Kaplan-Meier overall survival curve for THCApatients between high risk and low risk group;D:ROC curve showed the predictive value of risk signature.E,G,I:The IH Cexpression of IGF2BP1、IGF2BP2、RBM15 innormaltissue;F,H,J:The IH Cexpression of IGF2BP1、IGF2BP2、RBM15intumortissue.

有研究表明，在甲狀腺癌細胞，METTL3表達升高，可促進THCA的發(fā)生發(fā)展，但其未在組織水平中驗證[31]。本研究分析了m6ARNA甲基化調(diào)節(jié)基因THCA和正常組織中的表達，發(fā)現(xiàn)IGF2BP2、IGF2BP3、HNRNPC表達上調(diào)，而RBM15、METTL14、YTHDC1、FTO、ZC3H13、IGF2BP1、YTHDF3、YTHDC2、WTAP、HNRNPA2B1、 KIAA1429、 METTL3、 ALKBH5、YTHDF1表達下調(diào)。METTL3的表達水平與上述細胞學研究水平有所區(qū)別。表達上調(diào)的基因中，IGF2BP2和IGF2BP3均有文獻報道。IGF2BP2促進導致結(jié)直腸癌的發(fā)病機理和發(fā)展[32]。在胰腺癌中，IGF2BP3表達較高，且與患者的整體生存有關(guān)[33]。在卵巢透明細胞癌中，過表達的IGF2BP3是一個潛在的致癌基因[34]。肺腺癌中，HNRNPC表達上升[35]。表達下調(diào)的基因中，部分基因的文獻報道與本文一致。ALKBH5在結(jié)腸癌中被下調(diào)[36]。YTHDC2在直腸癌中是保護基因[37]。METTL3和YTHDF1與肺腺癌更好的OS和RFS相關(guān)[35]。過量表達METTL14在體外抑制HepG2細胞的遷移和侵襲，并在體內(nèi)抑制細胞增殖和轉(zhuǎn)移[38]。IGF2BP1基因表達情況在不同文獻報道自相矛盾。在HCC中，IGF2BP1抑制細胞增殖并促進細胞凋亡[39]。而在卵巢癌中，IGF2BP1升高且預后不良[40]。KIAA1429、RBM15 、HNRNPA2B1、FTO的文獻報道結(jié)果與本文相反。KIAA1429、RBM15 、HNRNPA2B1在肺腺癌中高表達[35]。同樣，F(xiàn)TO在乳腺癌中高水平表達，并與其不良預后相關(guān)[20]。究其原因，這可能與腫瘤異質(zhì)性有關(guān)，同一基因?qū)Σ煌┌Y的影響不同。

表1 風險分組臨床特征(n,%)Tab.1 The clinical in formation in different risk group

圖4 風險評分和臨床特征對THCA患者預后的影響Fig.4 Effects of the risk signature and clinical information on the prognosis of CESC patients.A:The heat map showed the different expression of m6A RNA methylation regulators and the distribution of clinical information between high risk group and low risk group;Univariate Cox regression(B)and Multivariate Cox regression (C) analysis of the correlation between clinical information and overall survival.

圖5 risk Score在THCA患者臨床亞組中的預測價值Fig.5 Prognostic value of the risk signature in THCA patients stratified into different clinical subgroup.Kaplan-Meiersurvival for patientin(A)N0and(B)stage3/4.

圖6 在驗證集中驗證預后模型Fig.6 Validation of the prognostic signature in tested group.A:Kaplan-Meier overall survival curve for THCA patients between high risk and low risk group in tested group;B:ROC curve showed the predictive value of risk signature in tested group;C:Univariate Cox regression analysis the correlation between clinical information and overall survival in tested group.*P＜0.05.

m6ARNA甲基化調(diào)節(jié)基因?qū)Π┌Y的早期診斷和治療產(chǎn)生重要影響，與癌癥的預后值也緊密聯(lián)系[41]。經(jīng)Lasso cox回歸分析，有4個m6ARNA甲基化調(diào)節(jié)基因IGF2BP2、FTO、IGF2BP1、RBM15，被篩選出以構(gòu)建風險模型。利用ROC分析預測風險評分，AUC為0.713，說明本研究的風險模型構(gòu)建良好。隨后按風險評分分為高危組和低風險組的CESC患者，經(jīng)單變量和多變量cox回歸分析顯示，風險評分是一個獨立的甲狀腺癌預測模型因素。此外，我們進一步評價了臨床特征因素中的風險評分，結(jié)果發(fā)現(xiàn)高風險和低風險組的OS在N0和stage3/4亞組中有顯著性差異，且高風險組OS均明顯更低。當然本研究有一定的局限性，即只分析了TCGA中的THCA數(shù)據(jù)，沒有進行臨床實驗驗證，這將是后期研究重點。

綜上所述，m6ARNA甲基化調(diào)節(jié)基因與THCA臨床特征密切相關(guān)。本研究選擇了4個m6ARNA甲基化基因，構(gòu)建了一個可以作為獨立預后因子的風險模型，為THCA中m6ARNA甲基化的后續(xù)研究提供了重要證據(jù)。