喉鱗狀細胞癌m6A相關基因預后預測模型的構建及分析

韓弈垣 王雪梅 賀 晴 曹曉林

喉癌是常見的頭頸部惡性腫瘤，喉鱗狀細胞癌(laryngeal squamous cell carcinoma，LSCC)是其最常見的病理類型，約占總發病數的95%[1]。我國喉癌發病數和死亡數占癌癥總發病和總死亡的比例分別是 0.67%和0.57%，標化發生率和病死率分別為1.3/10萬和0.7/10萬，且男性發生率和病死率遠高于女性[2,3]。在過去的10年中，盡管LSCC的診斷技術及治療方案不斷改進，其5年總生存率仍從66%下降至63%[4]。因此，尋找影響LSCC患者預后的生物靶點和分子機制越來越受到重視。

N6-甲基腺苷(N6-methyladenosine，m6A)是在真核生物中mRNA和lncRNA甲基化修飾的常見形式，受甲基轉移酶、脫甲基酶以及結合蛋白的動態調節，并參與許多重要的生理過程，例如干細胞分化和多能性，晝夜節律，胚胎發生和DNA損傷反應等[5]。目前已知的m6A相關基因包括YTHDC1、YTHDF2、YTHDC2、YTHDF1、FTO、ALKBH5、HNRNPC、METTL14、ZC3H13、KIAA1429、RBM15、WTAP和METTL3[6]。近年來有研究報道，異常的m6A修飾與頭頸部鱗狀細胞癌的發生及預后密切相關[7,8]。因此，本研究旨在探尋m6A相關基因在LSCC組織中的差異表達情況，并在此基礎上構建出LSCC的預后預測模型。

材料與方法

1.數據下載：從TCGA數據庫下載RNA-seq數據及相關臨床數據，篩選標準如下：①腫瘤原發部位：喉；②病理類型：鱗狀細胞癌；③腫瘤數據來源：TCGA-頭頸部鱗狀細胞癌；④數據類型：RNA-seq；⑤具備完整的臨床信息，包括生存時間、生存狀態、年齡、性別、腫瘤分級、臨床分期、腫瘤范圍、淋巴結轉移、遠處轉移。經篩選共納入123例RNA-seq數據樣本(包括111例腫瘤樣本及12例正常樣本)及92例具有完整臨床資料的LSCC患者臨床數據樣本。

2.m6A相關基因差異表達分析：利用Perl10.0.0軟件及ensembl數據庫將樣本中的ensembl id轉換為基因名，并整合各樣本RNA-seq數據，提取所有樣本中m6A相關基因的表達量。利用R4.0.3軟件中的Wilcoxon秩和檢驗計算腫瘤組織及正常組織中m6A相關基因(包括YTHDC1、YTHDF2、YTHDC2、YTHDF1、FTO、ALKBH5、HNRNPC、METTL14、ZC3H13、KIAA1429、RBM15、WTAP和METTL3)差異表達情況，并分別用pheatmap包、vioplot包進行熱圖及小提琴圖的繪制。以P<0.05為條件篩選LSCC樣本中顯著差異表達的m6A相關基因。

3.預后預測模型的構建及評估：為進一步進行變量選擇，避免m6A相關基因過度擬合，通過R4.0.3軟件中的survival包及glmnet包對顯著差異表達的m6A相關基因進行Lasso回歸分析及交叉驗證，并構建風險評分模型。根據該模型計算臨床樣本中各組患者的風險評分，并以風險評分中位值將患者分為高風險組和低風險組。利用survival包繪制Kaplan-Meier生存曲線，比較兩組間LSCC患者5年生存差異并進行Log-Rank檢驗。然后用survivalROC包繪制時間依賴的受試者工作特征曲線(receiver operating characteristic，ROC)評估此模型的預測性能。

4.獨立預后分析：通過R4.0.3軟件中的survival包及forestplot包對該預后預測模型進行單因素和多因素COX回歸分析，并繪制森林圖，從而評估各臨床變量(包括年齡、性別、腫瘤分級、臨床分期、T分期、N分期和M分期)及該模型的風險評分與患者預后的關系，進而判斷該模型能否作為獨立的預后指標。

結 果

1.m6A相關基因差異表達分析：LSCC組織樣本和正常組織比較后，共篩選出9個顯著差異表達的m6A相關基因，包括METTL3、YTHDF1、YTHDC2、ALKBH5、HNRNPC、RBM15、WTAP、FTO和KIAA1429。除了YTHDC2在LSCC組織中表達下調以外，其余8個基因在LSCC組織中表達均上調(圖1、圖2)。

圖1 m6A相關基因差異表達的熱圖與正常組織比較，*P<0.01，**P=0.000

圖2 m6A相關基因差異表達的小提琴圖

2.預后預測模型的構建：對上述9個顯著差異的m6A相關基因進行Lasso回歸分析，并使用交叉驗證建立模型，避免m6A相關基因過度擬合，結果顯示，當變量數為9時，模型的偏差值最小(圖3、圖4)。該預后預測模型的風險評分公式如下：risk score=(KIAA1429表達值×0.0757)+(HNRNPC表達值×0.0486)-(FTO表達值×0.1308)-(WTAP表達值×0.0055)-(RBM15表達值×0.2427)-(ALKBH5表達值×0.0015)-(YTHDC2表達值×0.1660)-(YTHDF1表達值×0.0444)-(METTL3表達值×0.2598)。

圖3 Lasso回歸算法篩選變量動態過程圖

圖4 交叉驗證參數λ選擇過程圖

3.預后預測模型的評估：整理從TCGA數據庫下載的臨床樣本數據，篩除信息不全的臨床樣本，共得到92例具有完整臨床資料的LSCC患者臨床數據樣本(表1)。根據風險評分公式計算92例LSSC患者的風險評分，并以風險評分中位值將患者分為高風險組和低風險組，其中高風險組36例，高風險組56例。Kaplan-Meier生存曲線結果顯示，高風險組的總體生存率顯著低于低風險組，且兩組的5年生存率之間比較差異有統計學意義(P<0.01，圖5)。時間依賴的ROC曲線分析結果顯示該模型在預測LSCC患者預后生存期方面具有良好的效果(AUC=0.728，圖6)。

表1 92例LSCC患者臨床資料信息

圖5 高、低風險組患者Kaplan-Meier生存曲線

圖6 預后預測模型的ROC曲線及AUC值

4.獨立預后分析：單因素及多因素COX回歸分析表明：m6A相關基因預后預測模型的風險評分與LSCC患者的總體生存率顯著相關(P=0.000)，且可作為獨立的預后因素(P=0.000)。此外，性別和淋巴結轉移情況同樣與LSCC患者預后顯著相關(P<0.05)，其中性別也可作為獨立的預后因素(P=0.000，圖7、圖8)。

圖7 各臨床特征及風險評分的單因素COX回歸分析森林圖

圖8 各臨床特征及風險評分的多因素COX回歸分析森林圖

討 論

LSCC的發生、發展涉及復雜的分子生物學改變，包括多種原癌基因的激活及抑癌基因的抑制，其具體的分子機制尚未完全闡明。近年來越來越多的研究證實異常的m6A修飾與頭頸部鱗狀細胞癌的發生與發展有關，并可依此來預測患者的生存預后，指導后續個體化治療方案的制訂[8～10]。

為了探尋m6A相關基因在LSCC組織中的差異表達情況，筆者從TCGA數據庫下載了LSCC組織的基因表達情況及相關臨床數據，并用Perl10.0.0軟件及R4.0.3軟件對數據進行了預處理及差異表達分析。筆者發現大多數m6A相關基因在LSCC組織中異常表達，包括8個顯著上調的基因(METTL3、YTHDF1、ALKBH5、HNRNPC、RBM15、WTAP、FTO和KIAA1429)和1個顯著下調的基因(YTHDC2)。有研究表明，METTL3高表達可導致腫瘤抑制因子SOCS2的m6A修飾增加，誘導SOCS2降解，促進肝細胞肝癌的發生；而KIAA1429也同樣在肝癌組織中高表達，且與肝癌患者的預后相關[11,12]。YTHDF1可通過Wnt/β-catenin途徑參與腸道惡性腫瘤的發展；同時也有研究證實YTHDF1與頭頸部鱗狀細胞癌患者的鐵代謝相關，并可能成為靶向治療的生物學標志物[10,13]。ALKBH5在膠質母細胞瘤和乳腺癌中表達上調，且與預后密切相關[14,15]。HNRNPC可通過EMT途徑促進口腔鱗狀細胞癌的惡性表達[16]。RBM15可作用于c-Myc抑制巨核細胞增殖，參與急性巨核細胞白血病的發展；而FTO、WTAP的異常表達可促進急性髓系白血病的發生[9,17,18]。YTHDC2作為抑癌基因，在頭頸部鱗狀細胞癌中顯著下調，且與多種免疫細胞的免疫浸潤水平相關[19]。

近年來，基于m6A相關基因的預后預測模型被用來評估和預測各種腫瘤患者的預后，指導后續個體化治療方案的制定[20～22]。本研究通過對9個在LSCC組織中顯著差異表達的m6A相關基因進行Lasso回歸分析，構建了LSCC的m6A相關基因的預后預測模型。Kaplan-Meier生存曲線顯示高風險組5年生存率明顯低于低風險組，且時間依賴的ROC曲線結果顯示該模型在預測LSCC患者預后生存期方面具有良好的效果(AUC值=0.728)。后續研究通過單因素及多因素COX回歸分析證明該模型可作為獨立的預后因子(P=0.000)。此外，本研究還發現性別也可作為獨立的預后因子(P=0.000)，且男性患者的預后欠佳。這可能與男性更多地暴露于煙酒等危險因素及性激素的復雜作用有關[23]。

綜上所述，本研究構建了基于9個m6A相關基因預測LSCC患者預后的風險評分模型。該模型具有較好的敏感度和特異性，且其可以作為獨立的預后因子評估LSCC患者預后，對制定合理及有效的個體化治療方案具有重要的參考價值。但該模型的應用仍需要大量臨床研究數據支持和大規模、多中心的循證醫學證據加以驗證。