王文龍，沈聰，孫博韜，白寧，李新營，陳嘉欣，彭婀敏

（中南大學(xué)湘雅醫(yī)院1.甲狀腺外科3.國際醫(yī)療部，湖南長沙410008；2.中南大學(xué)湘雅醫(yī)學(xué)院，湖南長沙410013）

過去幾十年，甲狀腺癌的發(fā)病率在全世界范圍內(nèi)迅速上升[1-2]。據(jù)國家癌癥中心數(shù)據(jù)顯示，甲狀腺癌已躋升至女性惡性腫瘤的第3 位，且發(fā)病出現(xiàn)年輕化趨勢，已成為我國30 歲以下女性增長最快的惡性腫瘤[3]。甲狀腺乳頭狀癌是最常見的病理亞型，占85%以上[4-5]。手術(shù)是其主要治療方式，規(guī)范化手術(shù)治療后10年生存率達90%以上[6]，但仍有1%～9%的患者在初次就診時就已出現(xiàn)遠處轉(zhuǎn)移，嚴(yán)重影響了患者預(yù)后[7]。因此，闡述甲狀腺癌的發(fā)病機制對改善預(yù)后至關(guān)重要。

m6A 甲基化修飾是指RNA 的腺苷酸的第6 位N處發(fā)生動態(tài)可逆的甲基化過程，在調(diào)控mRNA 的轉(zhuǎn)錄、翻譯、穩(wěn)定、剪切、成熟等發(fā)揮了重要的生物學(xué)作用，其功能主要由甲基轉(zhuǎn)移酶（writer）、去甲基化酶（eraser）和識別蛋白（reader）決定[8-10]。越來越多的證據(jù)表明m6A 甲基化修飾在腫瘤的發(fā)生發(fā)展扮演了重要的角色[9,11-13]。METTL3 通過YTHDF2 依賴的途徑上調(diào)抑癌基因SOCS2 的m6A 甲基化水平，從而加速SOCS2 mRNA 的降解。METTL3 過表達促進肝癌細胞的增殖和遷移，相反，敲低METTL3 顯著抑制腫瘤生長和轉(zhuǎn)移[14]。研究[15]揭示了IGF2BP3 在結(jié)直腸癌中高表達，下調(diào)IGF2BP3 可以通過VEGF 和CCND1 的m6A 修飾來抑制血管生成和DNA 復(fù)制。以上研究證實m6A 甲基化調(diào)控因子高度參與癌癥發(fā)生發(fā)展，具有良好的潛在預(yù)后價值。然而，目前對于m6A 甲基化調(diào)控因子在甲狀腺癌致病調(diào)控機制尚不清楚。

本研究利用癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）數(shù)據(jù)庫中甲狀腺癌RNA-sep數(shù)據(jù)，系統(tǒng)性分析了20 個m6A 甲基化調(diào)控因子在甲狀腺癌中的表達水平，再利用聚類分析和Lasso Cox 回歸分析篩選出預(yù)后相關(guān)基因及構(gòu)建風(fēng)險預(yù)測模型并進行危險分層，為甲狀腺癌患者的臨床診療和預(yù)后改善提供了重要的理論參考依據(jù)。

1 資料與方法

1.1 數(shù)據(jù)的獲取與處理

從TCGA 數(shù)據(jù)庫（https://portal.gdc.cancer.gov）下載了包含510 例甲狀腺癌組織和58 例正常組織樣本的RNA-seq 數(shù)據(jù)集及相應(yīng)臨床數(shù)據(jù)。臨床病理信息包括性別、年齡、腫瘤學(xué)多灶性、TNM 分期、淋巴結(jié)轉(zhuǎn)移、組織學(xué)類型、腫瘤位置和生存時間等，剔除RNA 序列以及臨床病理數(shù)據(jù)不全的樣本，共492 例癌組織和58 例正常組織樣本納入最終分析。所有患者的病理分期按照美國癌癥聯(lián)合委員會（AJCC）2018年1月頒布了第8 版甲狀腺癌TNM分期。根據(jù)之前m6A 甲基化在腫瘤中的研究[11,16]，共篩選20 個m6A 甲基化調(diào)控因子，包括甲基轉(zhuǎn)移酶（WTAP、METTL3、METTL14、METTL16、KIAA1429、ZC3H13 和 RBM15）、識 別 蛋 白（HNRNPA2B1、IGF2BP1/2/3、YTHDC1/2、FMR1、LRPPRC 和YTHDF1/2/3）和去甲基化酶（ALKBH5和FTO）。應(yīng)用R 語言的“Limma”軟件包分析上述基因在癌與正常樣本的差異表達。

1.2 聚類分析

為了從功能上闡明甲狀腺癌中m6A 調(diào)控因子的生物學(xué)特性，使用“consensusclusterplus”軟件包對甲狀腺癌患者進行一致性聚類分析并將患者分成cluster 1 和cluster 2 組。主成分分析（PCA）檢測基于m6A 調(diào)控因子表達是否能將兩組患者明顯區(qū)分，并探討不同的臨床病理因素和總體生存率在兩組聚類中是否存在差異。

1.3 風(fēng)險預(yù)測模型的建立與評估

為了探討m6A 調(diào)控因子在甲狀腺癌預(yù)后的作用，通過“glmnet”軟件包進行Lasso Cox 回歸分析篩選預(yù)后相關(guān)的基因，并計算基因的相關(guān)回歸系數(shù)，根據(jù)回歸系數(shù)加權(quán)建立甲狀腺癌預(yù)后風(fēng)險評估公式，risk score=Coei*Expi,Coei表示回歸系數(shù)，Expi表示m6A 調(diào)控因子的表達量。據(jù)此公式計算每個患者的風(fēng)險評分，以風(fēng)險評分的中位數(shù)將患者分為高危組和低危組，比較兩組總體生存率是否存在差異，ROC 曲線下面積（AUC）評估模型的預(yù)測能力。

1.4 統(tǒng)計學(xué)處理

采用R 3.6.0 軟件進行統(tǒng)計學(xué)分析，Wilcoxon 秩和檢驗分析20 個m6A 甲基化調(diào)控因子在腫瘤和正常組織的差異表達。計量資料用t檢驗，計數(shù)資料用χ2檢驗。采用Spearman 相關(guān)分析基因間的相關(guān)性，0.1≤|r|≤1.0 時定義為存在相關(guān)；生存曲線采用Kaplan-Meier 法和Log-rank 檢驗，P＜0.05 為差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 一般資料及m6A 甲基化調(diào)控因子的基因表達與相關(guān)性

本研究共納入492 例甲狀腺癌患者，平均年齡47.29 歲，男女比列1∶2.76；臨床病理TNM 分期：I 期276 例，II 期52 例；III 期109 例，IV 期55 例；T 分期：T1期135例，T2期162例，T3期172例，T4期23例；44.9%的患者出現(xiàn)頸部淋巴結(jié)轉(zhuǎn)移；9 例患者出現(xiàn)遠處轉(zhuǎn)移；228 例腫瘤呈現(xiàn)出多灶性，占46.3%；腫瘤位于峽部、左側(cè)、右側(cè)、雙側(cè)分別為22、172、214，84 例（表1）。

表1 492例甲狀腺癌患者基本臨床病理資料Table 1 The clinicopathologic information of the 492 patients with thyroid cancer

19 個m6A 甲基化調(diào)控因子在甲狀腺癌和正常組織表達具有統(tǒng)計學(xué)差異（P＜0.05），其中HNRNPC（P＜0.001)、IGF2BP2（P＜0.001）、FMR1（P＜0.05）在甲狀腺癌組織中明顯高表達。相反，ALKBH5（P＜0.001）、FTO（P＜0.001）、METTL3（P＜0.001）、METTL14（P＜0.001）、METTL16（P＜0.05）、WTAP（P＜0.001）、YTHDF1（P＜0.001）、YTHDF3（P＜0.001）、YTHDC2（P＜0.001）、YTHDC1（P＜0.001）、ZC3H13（P＜0.001）、HNRNPA2B1（P＜0.001）、RBM15（P＜0.001）、IGF2BP1（P＜0.001）、IGF2BP3（P＜0.001）和LRPPRC（P＜0.001）在甲狀腺癌標(biāo)本中表達明顯下調(diào)（圖1A）。此外，大部分m6A 甲基化調(diào)控因子都顯示正相關(guān)，其中，METTL14 與YTHDC1、YTHDF3與ZC3H13 呈現(xiàn)最強正相關(guān)（r=0.73）；而ALKBH5 與IGF2BP2呈明顯負相關(guān)（r=-0.33）（圖1B）。為了進一步了解19 個m6A 甲基化調(diào)控因子相互關(guān)系，PPI 構(gòu)建了網(wǎng)絡(luò)相互作用網(wǎng)（圖1C）。

圖1 m6A甲基化調(diào)控因子的基因表達及相關(guān)性A：20個m6A甲基化調(diào)控因子在甲狀腺癌和正常組織表達水平；B：Spearman相關(guān)分析m6A甲基化調(diào)控因子相關(guān)性；C：PPI網(wǎng)絡(luò)評估m(xù)6A甲基化調(diào)控因子相互作用關(guān)系Figure 1 Gene expression and correlations of the m6A methylation regulatorsA:The mRNA levels of the 20 m6A regulators in thyroid cancer; B:Spearman correlation analysis of co-expressions among the m6A regulators; C:PPI network of the m6A regulators

2.2 一致性聚類分析m6A甲基化調(diào)控因子

為了確定最佳的聚類數(shù)，采用一致性聚類分析對492 例甲狀腺癌患者進行聚類分組，共識矩陣k 值為2 時，甲狀腺癌樣本之間不存在交叉，聚類穩(wěn)定性好（圖2A）。因此，492 例甲狀腺癌患者被分成cluster 1 和cluster 2 組，主成分分析可以將兩組患者明顯區(qū)分（圖2B）。聚類熱圖示，cluster 1組頸淋巴結(jié)轉(zhuǎn)移發(fā)生率明顯高于cluster 2（P＜0.01），而TNM 分期、T 分期、M 分期、腫瘤多灶性、腫瘤位置、年齡、性別在兩組之間無差異（均P＞0.05）（圖2C）。生存分析示，cluster 1 生存時間低于cluster 2（P=0.048）（圖2D）。

圖2 一致性聚類分析m6A甲基化調(diào)控因子A：k值為2時，甲狀腺癌樣本之間不存在交叉，聚類穩(wěn)定性好；B：主成分分析區(qū)分cluster 1和cluster 2；C：不同的聚類與臨床病理的關(guān)系；D：總體生存率在兩組聚類的情況Figure 2 Consensus clustering analysis of the m6A RNA methylation regulatorsA:The most appropriate selection with clustering stability when k=2 used; B:Principal component analysis to distinguish cluster 1 and cluster 2; C:Heatmap and clinicopathologic features of the two clusters;D:The overall survival of cluster 1 and cluster 2

2.3 m6A甲基化調(diào)控因子構(gòu)建風(fēng)險預(yù)后模型及檢驗

Lasso Cox 回歸分析篩選4 個基因（IGF2BP2、RBM15、YTHDF1、YTHDF3），風(fēng)險評分=（-0.129×IGF2BP2）+（0.308×RBM15）+（0.331×YTHDF1）+（0.813×YTHDF3）（圖3A），根據(jù)回歸系數(shù)計算每個患者的風(fēng)險評分，以風(fēng)險評分的中位數(shù)將患者分為高危組和低危組，風(fēng)險評分和生存狀態(tài)分布如圖3B 所示，高風(fēng)險患者中死亡人數(shù)明顯增加。相比低風(fēng)險組患者，高風(fēng)險患者生存期明顯縮短（P=0.007）（圖3C），ROC 曲線下面積示該模型可預(yù)測甲狀腺癌患者的預(yù)后（AUC=0.731）（圖3D）。

圖3 m6A甲基化調(diào)控因子構(gòu)建風(fēng)險預(yù)后模型及檢驗A：Lasso Cox回歸模型的m6A基因及相關(guān)系數(shù)；B：風(fēng)險評分和生存狀態(tài)分布；C：高、低危組的生存曲線；D：ROC曲線評估模型的預(yù)測效能Figure 3 Construction od prognostic risk model and validationA:Coefficients calculated by Lasso Cox regression analysis;B:Risk score and survival status distribution;C:Kaplan-Meier analysis of overall survival in high-risk and low-risk groups;D:Assessment of the prediction accuracy by ROC curve

3 討論

甲狀腺癌是一種常見的內(nèi)分泌惡性腫瘤，發(fā)病率在世界范圍內(nèi)呈快速上升趨勢，絕大部分甲狀腺癌預(yù)后良好[17-18]。然而一些分化較差的甲狀腺癌和高危型甲狀腺乳頭狀癌具有高度侵襲性，易出現(xiàn)局部侵犯和遠處轉(zhuǎn)移，且尚無特效治療方法[19-20]。因此闡述甲狀腺癌發(fā)生和演變進展的分子機制對疾病的早期干預(yù)和預(yù)后預(yù)測具有重要的臨床應(yīng)用價值。

m6A 修飾是近年來的研究熱點，其作用主要參與轉(zhuǎn)錄后修飾，同時也參與了DNA 修復(fù)、細胞分化、細胞重編程、細胞應(yīng)激等生命活動[21-22]。m6A 甲基化在不同的腫瘤中均有報道，在甲狀腺癌中m6A 相關(guān)基因的作用仍不清楚，Ye 等[23]發(fā)現(xiàn)IGF2BP2 在甲狀腺癌中高表達，lncRNA MALAT1 通過競爭性結(jié)合miR-204，通過m6A 修飾識別，上調(diào)IGF2BP2，增強MYC 表達，從而刺激甲狀腺癌細胞增殖、遷移和侵襲，同時減弱腫瘤生長和細胞凋亡。然而該研究樣本量有限無法驗證lncRNA MALAT1/mir-204/IGF2BP2/m6A-MYC 軸在甲狀腺癌進展機制。Wang 等[24]報道METTL3 作為致癌基因通過上調(diào)TCF1 的m6A 甲基化水平促進甲狀腺癌的增值和侵襲。然而He 等[25]研究揭示METTL3 作為抑癌基因結(jié)合YTHDF2 通過c-Rel 和RelA 失活NF-κB 通路改變腫瘤相關(guān)的中性粒細胞浸潤來調(diào)節(jié)腫瘤生長。同樣，Hou 等[26]通過生物信息學(xué)證實METTL3 在甲狀腺癌組織中明顯低表達。上述研究證實m6A 相關(guān)基因表達水平和功能機制在甲狀腺癌中仍需進一步闡述。本研究初步驗證了METTL3 在正常組織中表達量明顯高于甲狀腺癌組織（P＜0.001），提示METTL3 可能是一種抑癌基因，可改善甲狀腺癌患者的預(yù)后，這也進一步驗證了He 等[25-26]的研究結(jié)果。

為了從功能上闡明甲狀腺癌中m6A 調(diào)控因子的生物學(xué)特性，本研究對492 例甲狀腺癌患者進行一致性聚類分析后并分成cluster 1 和cluster 2，cluster 1 生存期低于cluster 2（P＜0.05），相反，頸淋巴結(jié)轉(zhuǎn)移發(fā)生率明顯高于cluster 2（P＜0.01）。暗示m6A 調(diào)控因子能預(yù)測頸淋巴結(jié)轉(zhuǎn)移并影響甲狀腺癌患者預(yù)后，具有非常重要的臨床實際意義。聚類分析在肝癌[27]、腎癌[16]、乳腺癌[11]、胃癌[28]等惡性腫瘤預(yù)后分析中均起到較好的作用。

此外，為了探討m6A 調(diào)控因子在甲狀腺癌預(yù)后的作用，應(yīng)用Lasso Cox 回歸分析篩選的4 個基因（IGF2BP2、RBM15、YTHDF1、YTHDF3）構(gòu)建風(fēng)險評估模型，相比IGF2BP2，其余3 個基因在甲狀腺癌的研究中甚少，RBM15 不具有甲基轉(zhuǎn)移酶的催化功能，但它可以結(jié)合METTL3 和WTAP，并將這兩個蛋白引導(dǎo)到特定的RNA 位點進行m6A 修飾[8,29]；而YTHDF1 和YTHDF3 均屬于YTH 域家族，YTHDF1 通過與起始因子相互作用增強mRNA 翻譯和蛋白質(zhì)合成，YTHDF3 分別通過與YTHDF1 相互作用增強RNA 翻譯，與YTHDF2 結(jié)合促進RNA 降解[9,30]。相比低風(fēng)險組患者，高風(fēng)險患者生存期明顯縮短（P=0.007），ROC 曲線示該模型可以預(yù)測甲狀腺癌患者的預(yù)后（AUC=0.731）。相比TNM 分期和MACIS 預(yù)后評估系統(tǒng)，該模型基于RNA 測序數(shù)據(jù)，更符合當(dāng)今的精準(zhǔn)醫(yī)學(xué)理念。

同時，本研究也存在一些缺陷，首先，該研究基于生物信息學(xué)構(gòu)建的預(yù)后預(yù)測模型，無外部的實驗數(shù)據(jù)驗證，臨床實際應(yīng)用價值有限。其次，該研究只是初步探討了m6A 調(diào)控因子在甲狀腺癌預(yù)后作用，不能完整的闡述甲狀腺癌發(fā)生和演變的分子機制。最后，本研究基于492 例甲狀腺癌患者，樣本量不大，后續(xù)需要進行更大的獨立隊列研究來驗證。

綜上，本研究基于TCGA 數(shù)據(jù)庫分析甲狀腺癌m6A 調(diào)控因子表達水平，通過Lasso Cox 回歸分析篩選的4 個基因（IGF2BP2、RBM15、YTHDF1、YTHDF3）構(gòu)建風(fēng)險評估模型，AUC 提示該模型具有良好的預(yù)測效能，將為甲狀腺癌臨床診治提供潛在的理論依據(jù)。