基于TCGA數(shù)據(jù)庫構(gòu)建頭頸部鱗癌RNA甲基化調(diào)控因子預(yù)后模型

彭頤， 宗丹， 葛宜枝， 何俠

頭頸部鱗狀細(xì)胞癌(head and neck squamous cell carcinoma，HNSCC)起源于口腔、咽部和喉部的黏膜上皮，是全球第6大常見的惡性腫瘤。全球每年估計有890 000例新病例和450 000例死亡病例[1]。手術(shù)、化療、放療和免疫治療等多學(xué)科治療方案能夠顯著提高HNSCC患者的生存時間和生活質(zhì)量，但大多數(shù)患者會在3～5年內(nèi)迅速進(jìn)展，局部復(fù)發(fā)和遠(yuǎn)處轉(zhuǎn)移的風(fēng)險高達(dá)30%～60%，且復(fù)發(fā)和轉(zhuǎn)移后存活率較低[2]。因此，早期識別復(fù)發(fā)轉(zhuǎn)移高風(fēng)險人群，并給予積極的臨床干預(yù)有利于改善患者預(yù)后。

RNA甲基化是一種轉(zhuǎn)錄后水平調(diào)控基因表達(dá)的形式。其本質(zhì)是基因片段上的堿基位點被甲基化后，參與調(diào)控RNA的翻譯、剪切、表達(dá)與降解等生物學(xué)功能[3]。常見的RNA甲基化主要包括m6A、m5C、m1A與m7G等。RNA甲基化是一個可逆的過程，其修飾水平需依賴甲基轉(zhuǎn)移酶、結(jié)合蛋白和去甲基化酶等的動態(tài)調(diào)節(jié)[4]。近年來的研究表明，RNA甲基化在腫瘤的發(fā)展與轉(zhuǎn)歸中發(fā)揮著重要作用，不僅能通過多種細(xì)胞信號通路驅(qū)動細(xì)胞惡性轉(zhuǎn)化和化療耐藥等[4-5]，且與患者的預(yù)后密切相關(guān)[6]。因此，利用在腫瘤中異常表達(dá)的RNA甲基化調(diào)控因子，探索RNA甲基化與腫瘤惡性進(jìn)程之間的關(guān)聯(lián)，對于基礎(chǔ)研究或臨床治療腫瘤均有重要作用。

本研究旨在利用生物信息學(xué)尋找HNSCC組織中異常表達(dá)的RNA甲基化調(diào)控因子，探索這些基因?qū)︻A(yù)測患者預(yù)后的價值，并建立相應(yīng)的預(yù)后風(fēng)險模型，為臨床診斷及治療提供新思路。

1 資料與方法

1.1 數(shù)據(jù)收集 從TCGA數(shù)據(jù)庫(https://portal.gdc.cancer.gov/) 下載HNSCC患者腫瘤組織和正常頭頸組織的轉(zhuǎn)錄組數(shù)據(jù)以及患者臨床資料，包括年齡、性別、腫瘤分級與分期、生存時間與狀態(tài)等。應(yīng)用R統(tǒng)計軟件(R version 4.1.1)中的Limma包對測序所得的基因表達(dá)譜做歸一化處理。應(yīng)用caret包，將HNSCC數(shù)據(jù)隨機分成訓(xùn)練集與驗證集用于后續(xù)的數(shù)據(jù)處理。

1.2 差異表達(dá)的RNA甲基化調(diào)控因子的鑒定 通過查閱最新文獻(xiàn)，共選擇了45種已知的m6A、m5C、m1A 和m7G的 RNA 甲基化調(diào)控因子[7-9]。利用R軟件的Wilcoxon秩和檢驗鑒定在HNSCC樣本與正常樣本中差異表達(dá)的RNA甲基化調(diào)控因子(P<0.05，|log2 Fc|≥1)，用pheatmap與igraph包繪制表達(dá)量熱圖與差異基因相關(guān)性網(wǎng)絡(luò)圖。通過 STRING(https://string-db.org)構(gòu)建差異表達(dá)的RNA甲基化調(diào)控因子的蛋白相互作用(Protein-Protein Interaction Networks，PPI)網(wǎng)絡(luò)。

1.4 預(yù)后預(yù)測模型的評估 為了評估模型的預(yù)測性能，對構(gòu)建的模型進(jìn)行Kaplan-Meier生存分析、受試者工作特征曲線(receiver operating characteristic，ROC)、PCA主成分分析與t-SNE降維分析，并采用單因素及多因素COX回歸分析進(jìn)行預(yù)后分析。利用高、低風(fēng)險組間的差異表達(dá)基因，通過Cluster Profiler包進(jìn)行GO富集分析。最后，進(jìn)行了單樣本基因集富集分析(ssGSEA),比較高、低風(fēng)險組免疫細(xì)胞浸潤差異。

2 結(jié)果

2.1 RNA甲基化調(diào)控因子差異表達(dá)分析 從TCGA數(shù)據(jù)庫獲得44例正常頭頸組織以及502例HNSCC組織的轉(zhuǎn)錄組分析數(shù)據(jù)，并對45個RNA甲基化調(diào)控因子的表達(dá)水平進(jìn)行比較，共鑒定出38個差異表達(dá)的基因(圖1A)。其中，2個基因表達(dá)量在腫瘤組織中下降(YTHDC2、TET2)，36個基因呈現(xiàn)出上調(diào)趨勢。為進(jìn)一步探索差異表達(dá)的RNA甲基化調(diào)控因子的相互作用，我們進(jìn)行了蛋白質(zhì)-蛋白質(zhì)相互作用(PPI)網(wǎng)絡(luò)分析(圖1B)。設(shè)置PPI分析的最低交互評分為 0.7(高置信度)，并確定了排名前10的中心基因(METTL3、ALKBH5、VIRMA、METTL14、WTAP、FTO、RBM15、YTHDC1、ALKBH1、ZC3H13)。此外，本研究還分析了RNA甲基化調(diào)控因子的表達(dá)相關(guān)性，相關(guān)性閾值取0.4(中等程度相關(guān))(圖1C)。

1A：正常組織和腫瘤組織之間 RNA 甲基化調(diào)控因子表達(dá)水平的熱圖(藍(lán)色：低表達(dá)；紅色：高表達(dá)；*P<0.05， **P<0.01，***P<0.001)；1B：來自 STRING 的 RNA 甲基化調(diào)節(jié)因子的潛在相互作用；1C：RNA甲基化調(diào)控因子的共表達(dá)網(wǎng)絡(luò)

2.2 頭頸部鱗癌患者預(yù)后基因模型構(gòu)建 將HNSCC患者數(shù)據(jù)分為訓(xùn)練集與驗證集。在訓(xùn)練集中運用單因素COX回歸分析初步篩選與生存相關(guān)的RNA甲基化調(diào)控因子(圖2A)，確定了圖示的11個預(yù)后相關(guān)基因(P<0.2)。進(jìn)一步對上述11個差異有統(tǒng)計學(xué)意義的RNA甲基化調(diào)控因子進(jìn)行LASSO回歸分析(圖2B、2C)。根據(jù)λ值確定其中9個基因為模型的關(guān)鍵基因，并得出預(yù)后預(yù)測模型的風(fēng)險評分公式如下：風(fēng)險評分=[ExpYTHDC2×(-0.078)]+(ExpMETTL5×0.043)+(ExpHNRNPA2B1×0.016)+(ExpIGF2BP2×0.016)+(ExpIGF2BP1×0.031)+(ExpLRPPRC×0.005)+[ExpNSUN6×-(0.211)]+[ExpDNMT1×(-0.029)]+(ExpTRMT61B×0.107)。根據(jù)上述公式計算出的風(fēng)險評分，將患者平均分為低風(fēng)險組和高風(fēng)險組(圖2D)。主成分分析(PCA)和t-SNE分析顯示，不同風(fēng)險程度的兩組患者分布模式存在一定差別(圖2E、2F)。圖2G顯示了患者風(fēng)險評分與生存狀態(tài)的分布狀態(tài)(圖2G)。結(jié)果表明高風(fēng)險組患者相對于低風(fēng)險組死亡率更高，且生存時間更短。用Kaplan-Meier生存曲線評估該模型的預(yù)測能力，結(jié)果顯示風(fēng)險評分與HNSCC的生存相關(guān)，高風(fēng)險組的5年生存率僅有43%，而低風(fēng)險組的5年的生存率則相對較高為52%(P<0.001)(圖2H)。時間依賴的受試者工作特征曲線(ROC)顯示此風(fēng)險模型在預(yù)測生存狀態(tài)方面具有較好的診斷價值(1年AUC為0.619，3年AUC為0.653，5年AUC為0.601)(圖2I)。

2A：與OS相關(guān)的RNA 甲基化調(diào)控因子；2B：LASSO 回歸中參數(shù)選擇的交叉驗證；2C：LASSO 回歸篩選變量動態(tài)過程圖；2D：基于風(fēng)險評分的患者分布情況；2E：主成分分析；2F：t-SNE分析；2G：每位患者的生存狀態(tài)(低危人群：虛線左側(cè)；高危人群：虛線右側(cè))；2H：高風(fēng)險組和低風(fēng)險組患者 OS 的Kaplan-Meier 曲線；2I 預(yù)后模型的ROC曲線

2.3 預(yù)后模型的驗證 根據(jù)訓(xùn)練集建立的預(yù)后模型風(fēng)險評分公式計算驗證集中病例的風(fēng)險值，取中位風(fēng)險評分，將驗證集分為高風(fēng)險組與低風(fēng)險組(圖3A)。主成分分析(PCA)和t-SNE分析顯示高、低風(fēng)險組具有不同的分布模式(圖3B、3C)。低風(fēng)險組的患者比高風(fēng)險亞組中的患者具有更長的生存時間和更低的死亡率(圖3D)。Kaplan-Meier生存曲線還證明低風(fēng)險組的存活率明顯高于高風(fēng)險組的存活率(P<0.001)(圖3E)。ROC曲線分析顯示此模型有預(yù)測效力(1年AUC為0.614，3年AUC為0.649，5年AUC為0.673)(圖3F)。

3A：基于風(fēng)險評分的患者分布情況；3B：主成分分析；3C： t-SNE分析；3D：每位患者的生存狀態(tài)(低危人群：虛線左側(cè)；高危人群：虛線右側(cè))；3E：高風(fēng)險組和低風(fēng)險組患者 OS 的Kaplan-Meier 曲線；3F：預(yù)后模型的ROC曲線

2.4 風(fēng)險模式的獨立預(yù)后價值 單因素和多因素COX回歸分析來評估RNA甲基化調(diào)控因子預(yù)后預(yù)測模型。單因素COX回歸分析表明，該模型的風(fēng)險評分是預(yù)測TCGA隊列預(yù)后因素(HR=3.791，95%CI：1.477～9.726)(圖4A)。多因素分析同樣表明，在調(diào)整其他混雜因素后，風(fēng)險評分能夠作為HNSCC患者的獨立預(yù)后因素(HR=3.436，95%CI：1.299～9.093) (圖4B)。此外，患者的T分期、N分期與M分期情況同樣也可作為獨立的預(yù)后因素。本研究還做了高、低風(fēng)險組患者的臨床特征與預(yù)后預(yù)測模型中9個關(guān)鍵基因表達(dá)量的熱圖(圖4C)。結(jié)果顯示高、低風(fēng)險組的患者在腫瘤分期與T分期方面，差異有統(tǒng)計學(xué)意義(P<0.01)。

4A：TCGA 隊列風(fēng)險評分的單因素COX回歸分析；4B：TCGA 隊列風(fēng)險評分的多因素COX回歸分析；4C：臨床病理特征和風(fēng)險組之間聯(lián)系的熱圖(藍(lán)色：低表達(dá)；紅色：高表達(dá)；**P<0.01; ***P<0.001)

2.5 基于風(fēng)險模型的富集分析 在TCGA數(shù)據(jù)中，提取風(fēng)險模型分類中高、低風(fēng)險組之間差異表達(dá)的基因共2 630個。其中，1 576個基因在高危組中高表達(dá)，1 054個基因在高危組中低表達(dá)?；谶@些差異表達(dá)基因進(jìn)行GO富集分析。結(jié)果顯示，差異表達(dá)的基因主要參與到免疫應(yīng)答激活信號轉(zhuǎn)導(dǎo)、細(xì)胞吞噬作用以及免疫球蛋白相關(guān)的免疫應(yīng)答等機體免疫相關(guān)的生物學(xué)過程(圖5)。

注：條形越長表示富集的基因越多，紅色越深表示差異越明顯；q值表示調(diào)整后的P值

2.6 基于風(fēng)險模型的免疫功能比較 運用單樣本基因集富集分析(ssGSEA)，比較了HNSCC高、低風(fēng)險組免疫細(xì)胞的浸潤水平與免疫相關(guān)通路的差異。在TCGA數(shù)據(jù)中，高風(fēng)險組的免疫細(xì)胞浸潤水平普遍低于低風(fēng)險組(圖6A)。除巨噬細(xì)胞(Macrophages)浸潤水平降低不具有顯著性外，其余13種免疫細(xì)胞浸潤水平，差異有統(tǒng)計學(xué)意義(P<0.001)。免疫相關(guān)通路的差異分析也顯示，高風(fēng)險組相較于低風(fēng)險組存在普遍的免疫相關(guān)通路的抑制現(xiàn)象，有11種免疫相關(guān)通路均受到了不同程度的抑制(P<0.05)(圖6B)。

6A：低風(fēng)險組(藍(lán)框)和高風(fēng)險組(紅框)之間16種免疫細(xì)胞富集評分的比較；6B；低風(fēng)險組與高風(fēng)險組之間13種免疫相關(guān)通路類型富集評分的比較

3 討論

近年來，RNA甲基化在mRNA、tRNA、rRNA、snRNA等各種類型的RNA中的作用被廣泛關(guān)注，并日益成為腫瘤研究領(lǐng)域中的熱點。RNA甲基化水平的調(diào)控受到3類蛋白的影響：負(fù)責(zé)寫入甲基化的甲基轉(zhuǎn)移酶，能夠擦除甲基化的去甲基化酶，以及識別甲基化位點行使生物學(xué)功能的甲基化閱讀蛋白。多項研究證實RNA甲基化修飾能夠通過多種機制影響腫瘤的發(fā)生和發(fā)展，例如：參與腫瘤干細(xì)胞多功能性的形成與分化[10]，參與腫瘤細(xì)胞增殖[11]以及參與腫瘤細(xì)胞侵襲轉(zhuǎn)移[12]等，并且這些異常的RNA甲基化修飾的改變通常導(dǎo)致患者的不良預(yù)后[13]。此外，Li等[14]的研究證實，m6A甲基化修飾能夠通過靶向IL-7/STAT5/SOCS通路調(diào)控T細(xì)胞的穩(wěn)態(tài)，這意味著RNA甲基化在腫瘤免疫中也扮演著不可替代的角色。因此，本研究著重關(guān)注RNA甲基化修飾與腫瘤及腫瘤免疫的關(guān)系。

RNA甲基化是目前發(fā)現(xiàn)的170多種RNA轉(zhuǎn)錄后修飾中最主要的修飾方式之一，包括m6A、m5C、 m1A和m7G等。Chen等[15]的研究基于3個m6A調(diào)控因子(METTL4、METTL3、HNRNPC)建立了一個用于預(yù)測HNSCC預(yù)后的模型。該模型劃分的高、低風(fēng)險組能夠很好地估計患者的生存狀態(tài)，但是患者的風(fēng)險評分與其他臨床病理特征之間并沒有顯示出明顯的相關(guān)性。類似的文章還見于Han等[16]的研究，該研究基于m5C調(diào)控因子(NSUN5、DNMT1、DNMT3A)建立了HNSCC預(yù)后預(yù)測模型。然而，最近的一項研究表明在不同癌癥類型中，m6A與5mC調(diào)控因子之間存在表達(dá)的高度相關(guān)性，意味著m6A與5mC調(diào)控因子可能存在串?dāng)_現(xiàn)象。鑒于這種不同種類RNA甲基化修飾潛在的關(guān)聯(lián)，本研究選擇使用4種RNA甲基化的調(diào)控因子共同構(gòu)建預(yù)后模型。其中，RNA m6A甲基化是RNA最常見的轉(zhuǎn)錄后修飾，在本模型中也占據(jù)主導(dǎo)地位。m5C、m1A的調(diào)控因子也都參與到模型的構(gòu)建。而m7G的調(diào)控因子并未參與到風(fēng)險模型的構(gòu)建，考慮由于m7G甲基化的研究尚處于起步階段，目前已鑒定的調(diào)控因子過少導(dǎo)致。

本研究的預(yù)后模型涉及9個RNA甲基化調(diào)控因子，包括4個甲基轉(zhuǎn)移酶與5個結(jié)合蛋白。METTL5能夠催化18S rRNA的m6A甲基化并促進(jìn)腫瘤細(xì)胞生長與胚胎干細(xì)胞的分化[17-18]。TRMT61B則主要作為m1A甲基轉(zhuǎn)移酶催化線粒體的16S rRNA與tRNA[19]。DNMT1與NSUN6能夠催化m5C修飾，DNMT1已被證實在口咽癌、喉鱗癌等HNSCC過度表達(dá)，并且與不良預(yù)后相關(guān)[20-21]。NSUN6則是具有較強底物特異性的mRNA與tRNA m5C甲基轉(zhuǎn)移酶[22-23]。與本研究不同的是，NSUN6在胰腺癌中被鑒定為表達(dá)量下調(diào)，并可以通過調(diào)節(jié)細(xì)胞增殖抑制胰腺癌的發(fā)展[24]。NSUN6在HNSCC中的潛在生物學(xué)作用及機制還需進(jìn)一步研究來揭示。結(jié)合蛋白IGF2BPs家族在多項研究中被報道與mRNA的穩(wěn)定性和翻譯相關(guān)。IGF2BP1與IGF2BP2能夠分別通過介導(dǎo)BMI1的翻譯[25]與提高c-Myc的mRNA穩(wěn)定性[26]促進(jìn)口腔鱗癌的惡性進(jìn)展。hnRNP家族是另一類m6A結(jié)合蛋白，有報道認(rèn)為在HNSCC中HNRNPA2B1通過調(diào)節(jié)癌基因的剪接事件促進(jìn)上皮-間充質(zhì)轉(zhuǎn)化[27]。而高甲基化的LRPPRC在舌鱗狀細(xì)胞癌中被認(rèn)為具有一定診斷價值和預(yù)后預(yù)測潛力[28]。模型中比較特殊的m6A結(jié)合蛋白YTHDC2是為數(shù)不多的在HNSCC中表達(dá)量下降的調(diào)控因子，已被證實是與HNSCC預(yù)后、凋亡激活和泛素介導(dǎo)的蛋白水解有關(guān)的抑癌基因[29-30]。因此，本模型涉及的RNA甲基化調(diào)控因子與腫瘤的惡性進(jìn)展，尤其是HNSCC的增殖與侵襲存在著密不可分的關(guān)系，用此模型來預(yù)測HNSCC患者的預(yù)后是有據(jù)可循的。

有研究證據(jù)表明RNA甲基化修飾及其調(diào)控因子在病原體識別、免疫細(xì)胞激活和免疫反應(yīng)中起重要作用[31]。是免疫系統(tǒng)穩(wěn)態(tài)和激活的新型調(diào)節(jié)劑，也是腫瘤發(fā)生、轉(zhuǎn)移、治療抵抗和復(fù)發(fā)的切入點。本研究模型利用RNA甲基化調(diào)控因子表達(dá)量的差異劃分的高、低風(fēng)險組之間存在明顯的免疫功能與免疫浸潤差異，恰也證明了以上的觀點。

綜上所述，RNA甲基化調(diào)控因子在HNSCC中的異常表達(dá)可能在腫瘤的發(fā)生發(fā)展中起重要作用。本研究基于9個RNA甲基化相關(guān)調(diào)控因子構(gòu)建的預(yù)后模型，能有效對HNSCC患者進(jìn)行高、低風(fēng)險的分層，具有良好的預(yù)后預(yù)測效力，并且通過了驗證集的檢驗，能夠為臨床制定合理個體化治療方案提供參考。然而，本研究結(jié)果仍需要大量臨床數(shù)據(jù)的支持與多中心的循證醫(yī)學(xué)證據(jù)的驗證，RNA甲基化與機體免疫的作用相關(guān)性也有待進(jìn)一步發(fā)掘。