基于TCGA數據庫構建和評估口腔癌的免疫相關預后模型

巴 穎，張核子，余晨笛，盧曉萍，操利超

（深圳市核子基因科技有限公司，廣東 深圳 518071）

口腔癌（oral cancer）是頭頸部較常見的惡性腫瘤之一，已引起全球關注[1]。在國際分類定義中，口腔癌是包括嘴唇、舌頭、牙齦、口腔底部、頰粘膜、口腔前庭或磨牙后區等部位的癌癥[2]。盡管目前對口腔癌的診斷和治療方面已取得了一定的進展，但口腔癌患者的死亡率依然很高，主要原因之一是缺乏有效的預后生物標志物[3]。因此，探索口腔癌潛在的預后標志物和治療靶點是當前亟需解決的問題。免疫檢查點抑制劑（immune checkpoint inhibitors，ICIs）是一種很有前景的免疫療法，通過抑制負性調節受體發揮作用，如程序性細胞死亡受體1（programmed cell death protein 1，PD-1）和細胞毒性T淋巴細胞抗原4（cytotoxic T lymphocyte antigen 4，CTLA4），從而激活抗腫瘤免疫[4]。免疫治療在口腔癌中具有更好的療效，有關口腔癌免疫學的研究不斷增加[5，6]。然而，由于腫瘤免疫微環境的異質性和復雜性，只有一小部分患者受益于免疫治療，持續緩解的患者仍很少。許多免疫相關生物標志物被認為是口腔癌患者有用的預測指標，但目前尚缺乏一種理想的生物標志物用于臨床[7-10]。為了預測和改善口腔癌患者的預后，本研究基于癌癥基因組圖譜-頭頸部鱗狀細胞癌（the cancer genome atlas-head and neck squamous cell carcinoma，TCGA-HNSC）隊列構建和評估了口腔癌免疫相關預后模型，分析該預后模型與口腔癌患者臨床病理特征之間的關系，探索腫瘤免疫微環境的特征，包括腫瘤浸潤細胞構成、腫瘤突變負荷（tumor mutational burden，TMB）、以及PD-1/PD-L1/CLTA4的mRNA 表達水平，現報道如下。

1 材料與方法

1.1 數據獲取 癌癥基因組圖譜（the cancer genome atlas，TCGA）中mRNA 表達數據集和對應的臨床信息從UCSC Xena 平臺（https://xenabrowser.net/datapages/）下載，選擇隊列為癌癥基因組數據共享系統（genomic data commons，GDC）TCGA-HNSC，剔除非口腔癌相關部位的樣本，樣本信息見表1。從免疫學數據庫和分析平臺（the immunology database and analysis portal，ImmPort）（https://immport.niaid.nih.gov/）下載免疫相關基因列表，總計1509 個。

表1 TCGA數據集中樣本信息[n（%）]

1.2 構建免疫相關的預后風險模型 首先，利用R包對TCGA數據集進行差異基因分析，過濾標準為adjustedP-value＜0.05 和差異倍數大于1.5 倍（|Log2FC|＞0.585），然后與免疫基因列表取交集，得到免疫相關的差異基因。接著，使用R 包Survival 進行單變量Cox 比例風險回歸模型，選取與生存預后相關的差異基因。通過多因子回歸分析確定每個預后因子的回歸系數，建立預后風險評估模型，預測患者生存率。公式為：

風險分數=∑差異基因的回歸系數χi×歸一化處理后的基因表達量βi

1.3 繪制生存曲線和ROC 曲線 為了評估構建的風險預后模型的性能，利用R 包survivalROC 繪制受試者工作特征曲線（receiver operating characteristic curve，ROC），在ROC 曲線的轉折點選擇最佳風險評分臨界值，轉折點處真陽性和假陽性之間的差異最大。高于臨界值的患者屬于高風險評分組，低于臨界值的患者屬于低風險評分組，并使用R 包Survminer 繪制兩組的生存曲線。

1.4 構建和驗證列線圖 為了提高預后模型的性能，繪制列線圖，使用R 包rms，通過整合風險評分模型和臨床信息，包括年齡、性別和腫瘤分期，可視化不同患者特征的預后價值。此外，基于單因子回歸分析的森林圖分析臨床信息與總生存期（overall survival，OS）間的關系。其中，一致性指數（concordance index，C-index）可反映列線圖的預測準確性。

1.5 腫瘤免疫侵潤細胞類型的構成比較 基于構建的預后風險模型，將腫瘤樣本分為高風險評分組和低風評分組，并采用CIBERSORT 算法[11]，獲取22 種腫瘤浸潤免疫細胞的比例。隨后，通過非配對t檢驗比較高風險評分組和低風險評分組之間的免疫圖譜。

1.6 腫瘤免疫微環境的特征分析 使用R 包maftools可視化高風險評分組和低風險評分組的突變譜[12]。計算TMB 值，并通過非配對t檢驗統計高風險評分組和低風險評分組間TMB的差異；使用Kaplan-Meier 方法比較高風險評分組和低風險評分組間OS的差異；采用Wilcoxon 檢驗比較高風險評分組和低風險評分組間免疫檢查點及其配體的mRNA 表達水平。

2 結果

2.1 構建口腔癌的免疫相關預后模型 共得到1533個差異表達基因，其中上調差異基因數量為726 個，下調基因數量為807 個，見圖1A；與免疫基因列表取交集后，共得到73 個免疫相關基因，其中51 個基因下調，22 個基因上調。單因子回歸和多因子回歸分析表明，有6 個免疫相關的差異基因與OS 有關，見表2。其中，回歸系數見圖1B～圖1D。根據逐步回歸模型，Akaike 信息標準（Akaike information criterion，AIC）為1536.21，C-index 為0.63，見圖1E。

圖1 鑒定預后相關的差異表達基因

表2 與口腔癌預后相關的基因信息

2.2 預后模型的性能評估 基于已構建的預后風險模型，將口腔癌患者分為高風險評分組和低風險評分組，其中，cutoff 值設為0.02。可以看出，隨著風險得分的增加，生存時間呈現縮短的趨勢，并且高危組的死亡比例比低危組高，見圖2A、圖2B；基因CD40LG 和CMA1 在低風險評分組表達量高，在高風險評分組表達量低，而其余3 個基因趨勢相反，見圖2C；高風險評分患者的OS 比低評分患者預后較差，見圖2D。

圖2 預后風險評分分組和評估

2.3 預后模型的統計分析 繪制ROC 曲線和腫瘤分層分析進一步評估風險模型的性能，結果顯示，ROC曲線下面積（area under the curve，AUC）在3年時為0.678，4年時為0.671，5年時為0.683，見圖3A。Wilcoxon 檢驗表明，預后風險評分與N 分期（區域淋巴結轉移）密切相關（P=0.034），但與病理分期和T 分期無相關性（P＞0.05），見圖3B、圖3C、圖3D。

圖3 預后風險評估模型性能評價

2.4 構建和評估列線圖模型 在列線圖中，每個變量的得分可以在分數表上找到，然后通過計算總分來估計3、4 和5年的生存概率，見圖4A；森林圖顯示，患者年齡（＞60 歲）、腫瘤Ⅳ期和風險評分與OS 相關（P＜0.05），見圖4B。另外，繪制校準曲線來驗證列線圖的性能，結果顯示，構建的列線圖模型性能良好，預測曲線接近理想曲線，見圖4C～圖4E。此外，該列線圖（C-index：0.67）的預測準確性高于風險評分模型（C-index：0.63）。

圖4 列線圖模型的構建與驗證

2.5 體細胞突變統計及其與口腔癌生存預后的關系基于建立的免疫預后模型，將腫瘤樣本分為高風險評分組和低風險評分組，分別統計體細胞突變類型、突變模式和基因突變率。結果顯示，突變以錯義突變為主，主要突變類型為單核苷酸變異（single nucleotide variants，SNV），突變模式主要為C＞T。在高風險評分組中，突變率排名前10的基因為TP53、TTN、MUC16、FAT1、CDKN2A、KMT2D、LRP1B、NOTCH1、DNAH5 和SYNE1，在低風險評分組中，突變率排名前10的基因為TP53、TTN、FAT1、NOTCH1、CDKN2A、PIK3CA、CASP8、SYNE1、CSMD3 和LRP1B，見圖5A，圖5B；其中基因MUC16 和NOTCH1 在高風險評分組和低風險評分組中的突變率比較，差異有統計學意義（P＜0.05），見圖5C；對基因MUC16 和NOTCH1 分別進行生存曲線分析，兩個基因均不能作為獨立的預后因子（P＞0.05），見圖5D，圖5E；高風險評分組和低風險評分組的TMB 比較，差異無統計學意義（P＞0.05），見圖5F；但生存曲線分析顯示，TMB 值可作為獨立的預后因子，見圖5G。

圖5 體細胞突變分析統計及其與OS的關系

圖5 體細胞突變分析統計及其與OS的關系（續）

2.6 口腔癌腫瘤免疫微環境 基于CIBERSORT 算法，估算每個口腔癌患者中22 種免疫細胞的比例，并比較高風險評分組和低風險評分組之間的免疫細胞比例，結果顯示，有6 種免疫細胞類型存在差異（P＜0.05），見圖6A；另外，雖然免疫細胞NK.cells.activated 在高風險評分組和低風險評分組的免疫侵潤程度不存在差異，但其與生存預后相關（P＜0.05），見圖6B。Wilcoxon 檢驗結果表明，高風險評分組PD-L1（P=1.8×10-7）、PD-1（P=0.0079）和CTLA-4（P=9.7×10-7）的表達水平較低，見圖6C～圖6E。

圖6 口腔癌腫瘤免疫微環境

3 討論

免疫細胞在腫瘤的演變中起著重要作用[4]，目前ICIs 已應用于中晚期實體瘤的治療[13]。對于口腔癌，納武利尤單抗（nivolumab）和派姆單抗（pembrolizumab）均被批準用于治療復發/轉移性疾病[14]。然而，由于對腫瘤微環境特性的理解不足，這在一定程度上阻礙了免疫治療的廣泛應用。近年來，大量研究鑒定出了與口腔癌診斷和預后相關的免疫相關生物標志物[8，9，15，16]。但為了最大限度地發揮免疫治療的作用，還需要探索更多可靠的生物標志物。

本研究基于5 個免疫相關差異基因構建了口腔癌預后風險評估模型，這5 個基因為FGF9、CD40LG、CMA1、ESRRG、ADIPOQ。已有研究表明[7，17，18]，CD40LG、CMA1 和ESRRG 與口腔癌的預后顯著相關。FGF9 是一種生長因子，屬于FGFs 家族，FGF9是許多癌癥（包括肺腺癌、前列腺癌和肝癌）腫瘤演化的關鍵因子[19]。CD40LG 是一種跨膜蛋白，屬于腫瘤壞死因子超家族，存在于CD4+T 細胞、B 細胞、單核細胞和自然殺傷細胞及癌細胞中，具有增強腫瘤細胞的凋亡和抗原呈遞細胞的功能[17]。CMA1 是一種乳糜蛋白酶，由肥大細胞分泌，在癌癥中還沒有明確的研究。有研究表明[20]，通過增強ESRRG的表達和功能，可在癌細胞中發揮抑癌活性。ADIPOQ 表達的脂聯素可顯著抑制乳腺癌生長并誘導細胞凋亡[21]，但其在口腔癌中的分子機理還未見報道。

本研究結果顯示，免疫細胞NK.cells.activated的免疫侵潤程度和TMB 是口腔癌潛在的獨立預后標志。口腔癌高風險評分患者的免疫檢查點及其配體（PD-L1、PD-1 和CTLA-4）表達水平較低，這可能表明構建的免疫相關風險評分模型能夠為免疫治療提供一定的臨床指導意義。

綜上所述，本研究成功構建了免疫相關的口腔癌風險評估預后模型，并探索了腫瘤免疫微環境的特征，這可能有助于口腔癌患者的預后和免疫治療。