引用本文： 馮楊波, 熊延路, 趙晉波, 雷杰, 辛少偉, 喬天運, 周永生, 張曉, 姜濤, 韓勇. 全轉錄組分析建立食管癌風險預測模型及其評估. 中國胸心血管外科臨床雜志, 2023, 30(4): 576-585. doi: 10.7507/1007-4848.202102010 復制
全球食管癌（esophageal cancer，ESCA）的發病率在所有腫瘤中排第 7 位，總體死亡率居第 6 位，其 5 年生存率約 15%～25%[1-2]。在中國，ESCA 的發病率和死亡率在所有腫瘤中分別排第 5 位和第 4 位[3-4]。由于大多數 ESCA 診斷時已處于晚期，患者總體預后較差[2, 5]。2018 年，全球約有 572 000 例 ESCA 確診病例，509 000 人死于 ESCA。ESCA 的發病率和組織類型因地理位置的不同存在很大差異。ESCA 是一種異質性疾病，主要分為鱗狀細胞癌和腺癌兩種亞型。食管鱗狀細胞癌（esophageal squamous cell carcinoma，ESCC）是世界范圍內最常見的 ESCA 組織學亞型，主要見于東亞和中東地區，它的主要病因是大量飲酒和吸煙[6]。吸煙和飲酒產生的協同作用會顯著增加 ESCC 發生風險。此外，代謝酒精的酶突變、社會經濟水平偏低、口腔衛生差、營養缺乏也與 ESCC的發生有關[2]。食管腺癌（esophageal adenocarcinoma，EAC）主要見于高收入國家，Barrett’s 食管、超重和胃食管反流病是其主要危險因素[7-8]。
ESCA 的預后與腫瘤分期、腫瘤亞部位、組織學、患者的身體狀況及有無并發癥等相關[9-10]。目前，組織學分型和 TNM 分期仍是評估 ESCA 預后的主要方法。然而，現有指標無法充分展示患者全部疾病信息，也無法精確判斷患者預后。因此，建立有效的 ESCA 預后模型對完善 ESCA 預后評估系統、幫助臨床醫生優化治療方案及改善患者預后具有重大意義。
眾所周知，基因組的不穩定性和突變是腫瘤的重要特征，癌基因和抑癌基因的異常是 ESCA 發生發展的重要原因[11]。目前，對 ESCA 致癌基因的研究顯著推動了其治療的進步。首先，是針對驅動基因改變的分子靶向治療。例如，EAC 傾向于過表達血管內皮生長因子（vascular endothelial growth factor，VEGF）形成新生血管發揮促腫瘤進展作用。針對 VEGF 信號通路的靶向治療已經成為晚期 EAC 的二線治療方案[12]。此外，EAC 可出現人類表皮生長因子受體-2 增殖相關通路的活化突變。無論是作為單一療法還是與 Abraxane 聯用，ramucirumab 已被美國食品藥品監督管理局（Food and Drug Administration，FDA）批準用于 EAC 的二線治療[13]。其次，針對免疫檢查點基因的免疫治療也是 ESCA 治療的方向和熱點。例如，針對程序性死亡分子 1（programmed death 1，PD-1）的單克隆抗體 pembrolizumab 已被 FDA 批準為聯合陽性評分≥10%的復發性、局部晚期或轉移性程序性死亡分子 1 配體（PD-1 ligand，PD-L1）陽性的 ESCC 的二線選擇。隨著組學技術的飛速發展和生物信息學技術的廣泛應用，建立基于基因標簽的臨床預測模型已成為構建新的腫瘤預后評估系統的趨勢和潮流。
本研究從癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）中獲取了 ESCC 和 EAC 的 RNA 測序數據和相應的臨床數據，建立了基于關鍵基因的 ESCC 和 EAC 的風險預測模型。隨后，我們對模型的準確性進行了評估。最后，我們分析了模型的臨床意義及其與腫瘤免疫微環境的關系。這些工作將為準確評估 ESCC 和 EAC 預后、優化治療手段提供重要依據。
1 資料與方法
1.1 數據資料
從 TCGA 數據庫（
1.2 統計學分析
利用“edgeR”包進行差異分析來獲取 ESCA 的差異基因［adjusted P value<0.05，|log2（fold change）| > 1］。通過單變量 Cox 回歸分析篩選出與總生存期（overall survival，OS）相關的差異基因。利用“glmnet”包進行最小絕對收縮和選擇算子（least absolute shrinkage and selection operator，LASSO）回歸分析，進一步明確 OS 相關核心基因。通過“survival”、“survminer”包進行多變量 Cox 回歸分析，建立 ESCC 和 EAC 的風險評分預后模型。以風險評分中位數為截斷值，將 ESCC 和 EAC 患者分別分成高、低風險組。利用“pheatmap”包繪制參與風險評分基因的熱圖。繪制 Kaplan-Meier 曲線并進行 log-rank 檢驗，探討風險評分與患者生存狀態的關系。通過“survival ROC”包進行受試者工作特征（receiver operating characteristic，ROC）曲線分析，評估基于風險評分的生存預測的靈敏度和特異度。ROC 曲線下面積代表了預后準確性。對臨床特征進行二分類，基于該分類的風險評分差異分析通過 Wilcoxon 檢驗實現。利用“rms”和“survival”包，構建了含臨床分期、T 分期、N 分期、風險評分參數的列線圖。通過校準曲線評估實際生存率和預測生存率的一致性。計算 C-index，評估預后預測模型的準確性。利用“estimate”包估計每個樣本在腫瘤微環境中的免疫成分和基質成分含量。基因集富集分析（gene set enrichment analysis，GSEA）闡述高風險組和低風險組之間差異表達的基因所富集的生物學過程。從 MSigDB 數據庫下載 KEGG 作為 GSEA 分析的參考基因集，只有 NOM P<0.05 的富集通路才被認為有意義。利用 CIBERSORT 估計腫瘤樣本的腫瘤浸潤免疫細胞豐度，僅 P<0.05 的腫瘤樣品被選擇出來。高、低風險組之間免疫細胞差異和免疫檢查點基因表達的差異分析均采用 Wilcoxon 檢驗。P≤0.05 為差異有統計學意義。所有統計學分析均利用 R 語言完成。
2 結果
2.1 ESCC 和 EAC 生存相關核心基因的識別
腫瘤差異基因在腫瘤發生發展中發揮重要作用，對腫瘤差異基因的識別有助于加深我們對腫瘤的認識。我們從 TCGA 數據庫獲取了 ESCA 的 RNA 測序數據，利用“edgeR”包得到了 ESCA 的差異基因（圖 1a）。這些差異基因包括 1 482 個上調基因和 1 765 個下調基因。隨后，我們在 ESCC 和 EAC 中對 ESCA 差異基因的表達進行了單變量 Cox 回歸分析。在 ESCC 中，180 個基因和 OS 相關（P<0.05）；在 EAC 中，371 個基因和 OS 相關（P<0.05，圖 1b）。考慮到我們得到的相關基因比較多，為進一步篩選關鍵基因，我們進行了 LASSO 回歸分析并計算了回歸系數（圖 1c～d）。在 ESCC 中，我們找到 3 個與 OS 顯著相關的基因（NUDCD1、NUP155、HLA-B），在 EAC 中，我們找到 6 個與 OS 顯著相關的基因（CIRBP、SLC26A9、MSH2、HIST1H3F、PMAIP1、HIST1H3J）。
圖1
ESCA 差異基因的獲取及 OS 相關基因識別
a：TCGA 數據庫中 ESCA 差異基因；b：單變量 Cox 回歸分析確定 ESCC 和 EAC 的 OS 相關基因；c：通過 LASSO 回歸分析進一步確定與 ESCC 的 OS 相關的 3 個關鍵基因；d：通過 LASSO 回歸分析進一步確定與 EAC 的 OS 相關的 6 個關鍵基因；TCGA：癌癥基因組圖譜；ESCA：食管癌；ESCC：食管鱗狀細胞癌；EAC：食管腺癌；OS：總生存期；LASSO：最小絕對收縮和選擇算子；FC：差異倍數；DEGs：差異表達基因
2.2 ESCC 和 EAC 風險預測模型的構建和評估
我們利用前期獲得的生存相關核心基因表達值進行了多變量 Cox 回歸分析（逐步回歸，表 1～2），并依托“survival”包中的“predict（）”函數算法得到了每個樣本的風險評分值。其中，ESCC 得到了由 3 個差異基因構成的風險評分預測模型，EAC 得到了由 5 個差異基因構成的風險評分預測模型。這些風險評分模型主要用于預測 ESCC 和 EAC 患者的預后。
表1
食管鱗狀細胞癌風險預測模型所涉及的基因
表2
食管腺癌風險預測模型所涉及的基因
以風險評分的中位數為截斷值，將 ESCC 和 EAC 患者分為高、低風險組。無論是 ESCC 患者還是 EAC 患者，高風險組患者存活人數更少（圖 2a～b）。
圖2
風險評分預后模型的構建和驗證
a：風險評分、患者生存時間的分布，ESCC 患者生存狀態和 ESCC 患者預后相關基因表達譜的熱圖；b：風險評分、患者生存時間的分布，EAC 患者生存狀態和 EAC 患者預后相關基因表達譜的熱圖；c：以風險評分為分層因素的 ESCC 患者的 Kaplan-Meier 生存分析；d：ROC 曲線分析顯示風險評分在預測 ESCC 患者生存時間中的準確性；e：以風險評分為分層因素的 EAC 患者的 Kaplan-Meier 生存分析；f：ROC 曲線分析顯示風險評分在預測 EAC 患者生存時間中的準確性；TCGA：癌癥基因組圖譜；ESCC：食管鱗狀細胞癌；EAC：食管腺癌；ROC：受試者工作特征；AUC：曲線下面積
風險評分在預測 ESCC 和 EAC 患者預后上所體現的價值需要進一步評估。單變量 Cox 回歸分析表明，ESCC 和 EAC 的風險評分均與患者 OS 顯著相關［HR=1.221，95%CI（1.116，1.336），P<0.001；HR=1.055，95%CI（1.023，1.087），P<0.001］。我們納入風險評分、TNM 分期和臨床分期進行了多變量 Cox 回歸分析。結果顯示，ESCC 和 EAC 的風險評分與 OS 顯著相關［HR=1.318，95%CI（1.173，1.482），P<0.001；HR=1.048，95%CI（1.016，1.082），P=0.003］。這表明，風險評分可作為 ESCC 和 EAC 患者的獨立預后因子。隨后，我們采用 Kaplan-Meier 法和 log-rank 檢驗進行生存分析，評估風險評分水平與患者 OS 的關系。結果顯示，高風險的 ESCC 和 EAC 患者 OS 顯著短于低風險組（P=0.009，P<0.001）。為評估風險預測模型的預測效果，我們進行了時間依賴的 ROC 曲線分析。ESCC 和 EAC 的風險評分預測1年生存率的曲線下面積（area under the curve，AUC）分別為 0.765 和 0.810，這提示風險預測模型預測效果較好（圖 2c～f）。
2.3 ESCC 和 EAC 風險預測模型的臨床意義
為進一步闡明風險評分與臨床特征的關系，我們將臨床特征進行二分類，通過 Wilcoxon 檢驗比較了這些臨床特征的風險評分差異。在 ESCC 中，T 分期為 T1～2 的患者風險評分更高（P=0.039），而性別（P=0.606）、淋巴結轉移與否（P=0.457）、有無遠處轉移（P=0.224）、腫瘤分期的高低（P=0.480）與風險評分的關系皆不顯著。在 EAC 中，淋巴結轉移（P=0.013）、存在遠處轉移（P=0.031）、臨床分期為 Ⅲ～Ⅳ期（P=0.016）患者風險評分更高，而 T 分期（P=0.965）、性別（P=0.214）與風險評分無關。
為預測 ESCC 和 EAC 患者的 OS，我們利用“rms”和“survival”包構建了包括風險評分和臨床特征（臨床分期、T 分期、N 分期）的列線圖。ESCC 和 EAC 患者列線圖的 C-index 分別為 0.870（SE=0.050）和 0.807（SE=0.046），展現出了較高的辨識能力。隨后，利用校準曲線對實際生存期和預測生存期之間的一致性進行評估。校準曲線顯示出了較好的預測能力（圖 3）。
圖3
ESCC 和 EAC 風險評分的臨床意義
a：以風險評分、臨床分期、T 分期、N 分期參數構建列線圖預測 ESCC 患者的 1 年生存率；b：通過校準曲線評價 ESCC 患者的預測 1 年生存率和實際 1 年生存率的一致性；c：以風險評分、臨床分期、T 分期、N 分期參數構建列線圖預測 EAC 患者的 1 年生存率；d：通過校準曲線評價 EAC 患者的預測 1 年生存率和實際 1 年生存率的一致性；TCGA：癌癥基因組圖譜；ESCC：食管鱗狀細胞癌；EAC：食管腺癌
2.4 風險預測模型與腫瘤微環境的關系
腫瘤微環境的結構成分主要由基質細胞和募集的免疫細胞構成。我們利用“estimate”算法對 ESCC 和 EAC 樣本中免疫成分和基質成分的數量進行了估計。為揭示腫瘤樣本中基質成分比例和免疫成分比例與風險評分的關系，采用 Wilcoxon 檢驗進行差異分析。結果表明，ESCC 中免疫成分比例更高的樣本風險評分更高（P=0.001），而基質成分比例高、低組風險評分差異無統計學意義（P=0.133）。在 EAC 樣本中，基質成分（P=0.387）和免疫成分（P=0.051）比例高、低組的風險評分差異均無統計學意義。
隨后，我們分別對高、低風險組的 ESCC 和 EAC 樣本進行了 GSEA 分析。結果顯示，低風險組的 ESCC 顯著富集在抗原呈遞和自然殺傷細胞介導的細胞毒性生物學過程（圖 4a～b）。這意味著風險評分越高的腫瘤部分免疫功能受到抑制。然而，基于風險評分分組的 EAC 與免疫相關的生物學過程關系并不密切，這與我們前面的結果吻合。
為進一步揭示風險評分與腫瘤浸潤免疫細胞的關系，我們利用 CIBERSORT 算法分析了 ESCA 組織中腫瘤浸潤免疫亞群的比例，并在 ESCC 和 EAC 樣本中分別構建了 21 種、20 種免疫細胞的表達譜（圖 4c～d）。隨后，我們評估了高低風險組的免疫細胞表達是否存在差異。結果表明，ESCC 高、低風險組樣本中活化的肥大細胞差異有統計學意義（P=0.035），EAC 高、低風險組樣品中免疫細胞差異無統計學意義（P>0.05，圖 4e～f）。隨后，我們探討了基于高、低風險分組的常見免疫檢查點基因的表達特征。我們發現，ESCC 患者中，大部分免疫檢查點基因在高風險組中高表達（P<0.05，圖 4g）。然而，在 EAC 患者中高、低風險組的免疫檢查點基因表達差異無統計學意義（P>0.05，圖 4h）。
3 討論
ESCA 是一種致命性惡性腫瘤，給家庭和社會帶來了沉重的負擔。ESCA 預后判斷的準確性影響著治療方法的選擇和患者預后的改善。因此，建立一套準確可行的預后判斷體系具有重要意義。腫瘤細胞基因表達的異常和當前組學技術的發展為建立基于基因的預后判斷模型提供了理論和技術支撐[11, 14]。
本研究中，我們利用 ESCA 腫瘤樣本和正常樣本的 RNA 測序數據進行差異分析，得到 1 482 個上調基因和 1 765 個下調基因。隨后，我們進行單變量 Cox 回歸分析，在這些差異基因中找到了 180 個與 ESCC 的 OS 相關基因，找到了 371 個與 EAC 的 OS 相關基因。由于篩選得到基因較多，我們進行了 LASSO 回歸分析、多變量 Cox 回歸分析，找到 3 個與 ESCC 的 OS 顯著相關基因（NUDCD1、NUP155、HLA-B），6 個與 EAC 的 OS 顯著相關基因（CIRBP、SLC26A9、MSH2、HIST1H3F、PMAIP1、HIST1H3J），并依托“survival”包中“predict（）”函數得到 ESCC 和 EAC 的風險預測模型。基于該風險預測模型的分組準確地將患者 OS 的長短區分出來。針對 1 年生存率的預測，ESCC 和 EAC 風險預測模型的 AUC 分別為 0.765 和 0.810，顯示出了較好的預測效果。此外，該風險模型與臨床特征關系密切。ESCC 的風險預測模型與腫瘤 T 分期相關，EAC 的風險預測模型與淋巴結轉移、遠處轉移和臨床分期相關。隨后，我們構建了包含風險預測模型和相關臨床特征的列線圖。ESCC 的風險預測模型與腫瘤樣本的免疫成分相關，這與我們的 GSEA 分析結果相吻合。通過進一步分析 ESCC 風險預測模型和免疫細胞的關系，我們發現 ESCC 的風險預測模型和肥大細胞激活途徑相關。我們以風險預測模型作為分組依據，評估了免疫檢查點基因的表達情況。我們發現，大部分免疫檢查點基因的表達和 ESCC 風險預測模型顯著相關。這進一步證實了免疫逃逸在腫瘤進展中的重要作用。
事實上，該評分模型所涉及的基因具有重要的生物學意義，尤其在腫瘤的轉化和惡性進展中。如NUDCD1是一種在多種腫瘤中高表達或激活的被稱為腫瘤抗原的癌蛋白[15-16]。NUDCD1 參與腫瘤細胞的增殖、侵襲、凋亡和上皮間質轉化等生物學過程，發揮促癌作用[16-17]。NUP155在肝細胞癌中高表達，與 TP53 突變和細胞周期調節基因的過表達相關，其高表達患者往往預后較差[18-19]。一項關于 ESCC 的研究[20]表明，NUP155 在 ESCC 中高表達，但是對 NUP155 進行的生存分析結果卻顯示高表達者預后更好。NUP155 在 ESCC 中具體作用并未闡明，仍需要進一步研究。HLA-Ⅰ類分子（HLA-A、HLA-B、HLA-C）的高表達與 ESCA、結直腸癌和肝細胞癌的較好預后有關，這可能歸因于腫瘤對 CD8+T 細胞抗腫瘤免疫力比較敏感[21-23]。而在胃癌、肺癌中，HLA-Ⅰ類分子高表達預示著更差的預后，可能是因為腫瘤從自然殺傷細胞中逃逸，這一過程被稱為“自我識別缺失”[24-25]。Zhang 等[26]認為，HLA-Ⅰ類分子表達的下調可能是 ESCC 的獨立不良預后因素。當然，有的研究[27]也認為 HLA-Ⅰ類分子表達與 ESCC 預后無關。Jiang 等[28]認為，HLA-B 在 ESCC 中上調，并在 ESCC 中發現了1個新的融合基因（HLA-E 和 HLA-B）。而我們的研究顯示，HLA-B 在 ESCC 中高表達，與預后不良相關。這與以前的研究存在一定的區別，這可能與腫瘤異質性或者腫瘤標本免疫細胞浸潤程度有關。CIRBP 屬于富含甘氨酸的 RNA 結合蛋白家族，在應激條件時發揮關鍵作用[29]。一方面，CIRBP 可以發揮抑癌作用，如抑制腫瘤細胞增殖[30-31]。另一方面，CIRBP 可發揮促癌作用，如促進腫瘤細胞增殖、減少凋亡和衰老、參與上皮間質轉化、增強侵襲和遷移等過程[29,32-35]。我們的研究表明，CIRBP 在 ESCA 中低表達，其表達值越高預示著 EAC 預后越好。CIRBP 在 EAC 中的作用暫不明確，需要深入探究。SLC26A9 主要作為 Cl?/HCO3?交換器、Cl?通道和 Na+耦合的轉運蛋白而發揮作用[36-37]。基于 TCGA 數據的分析表明，SLC26A9 在肺鱗狀細胞癌中低表達，但在基底乳腺癌中高表達[38]。我們的研究表明，SLC26A9 在 ESCA 中低表達，扮演保護基因的角色。但是，SLC26A9 在 EAC 中的具體作用還需進一步研究。MSH2 是第1個被鑒定的錯配修復（mismatch repair，MMR）基因，其蛋白功能是識別錯配的堿基[39-40]。MSH2 表達的缺失在 Barrett 相關 EAC 中較少見，它是部分 Barrett 相關 EAC 的 MMR 缺陷原因之一[41]。Evans 等[42]的研究表明，78% 的 EAC 中存在 MSH2 表達。我們的研究認為 MSH2 在 ESCA 中上調，這與之前的研究并不矛盾。HIST1H3F屬于依賴復制的組蛋白基因家族。HIST1H3F 在侵襲性乳腺癌中高表達，其高表達與無復發生存期、無遠處轉移生存期較差有關[43]。同樣，我們的研究認為 HIST1H3F 在 ESCA 中高表達，與預后不良有關。PMAIP1 是 Bcl-2 蛋白家族的促凋亡成員，其表達可激活線粒體凋亡級聯反應并誘導氧化應激和促進鈣釋放[44]。然而，PMAIP1 常在黑色素瘤中高表達，環磷腺苷效應元件結合蛋白介導的 MEK/ERK 信號通路激活可促進 PMAIP1 的表達以促進自噬。PMAIP1 介導的自噬激活在營養缺乏條件下對黑色素瘤細胞有間接的抗凋亡作用[45]。PMAIP1 在 EAC 中的作用需要進一步探索。
當然，我們的研究也存在一些缺陷。首先，我們的研究基于轉錄組數據，轉錄水平并不完全等同于蛋白水平。其次，轉錄組數據主要包含蛋白編碼基因，而非編碼 RNA 在腫瘤中亦發揮重要作用[46]。最后，我們的研究缺乏詳細的機制研究。這些工作將在以后的研究中展開。
利益沖突：無。
作者貢獻：馮楊波負責醞釀和設計實驗，實施研究，采集、分析和解釋數據，起草文章；熊延路負責醞釀和設計實驗，研究經費支持，分析和解釋數據，對文章內容作批評性審閱；趙晉波負責對文章內容作批評性審閱，統計學分析，指導研究；雷杰負責對文章內容作批評性審閱；辛少偉、喬天運負責采集數據，對文章內容作批評性審閱；周永生負責分析和解釋數據，對文章內容作批評性審閱；張曉負責統計分析指導，對文章內容作批評性審閱；姜濤負責對文章內容作批評性審閱；韓勇負責醞釀和設計實驗，研究經費支持，指導研究，對文章內容作批評性審閱及定稿。
全球食管癌（esophageal cancer，ESCA）的發病率在所有腫瘤中排第 7 位，總體死亡率居第 6 位，其 5 年生存率約 15%～25%[1-2]。在中國，ESCA 的發病率和死亡率在所有腫瘤中分別排第 5 位和第 4 位[3-4]。由于大多數 ESCA 診斷時已處于晚期，患者總體預后較差[2, 5]。2018 年，全球約有 572 000 例 ESCA 確診病例，509 000 人死于 ESCA。ESCA 的發病率和組織類型因地理位置的不同存在很大差異。ESCA 是一種異質性疾病，主要分為鱗狀細胞癌和腺癌兩種亞型。食管鱗狀細胞癌（esophageal squamous cell carcinoma，ESCC）是世界范圍內最常見的 ESCA 組織學亞型，主要見于東亞和中東地區，它的主要病因是大量飲酒和吸煙[6]。吸煙和飲酒產生的協同作用會顯著增加 ESCC 發生風險。此外，代謝酒精的酶突變、社會經濟水平偏低、口腔衛生差、營養缺乏也與 ESCC的發生有關[2]。食管腺癌（esophageal adenocarcinoma，EAC）主要見于高收入國家，Barrett’s 食管、超重和胃食管反流病是其主要危險因素[7-8]。
ESCA 的預后與腫瘤分期、腫瘤亞部位、組織學、患者的身體狀況及有無并發癥等相關[9-10]。目前，組織學分型和 TNM 分期仍是評估 ESCA 預后的主要方法。然而，現有指標無法充分展示患者全部疾病信息，也無法精確判斷患者預后。因此，建立有效的 ESCA 預后模型對完善 ESCA 預后評估系統、幫助臨床醫生優化治療方案及改善患者預后具有重大意義。
眾所周知，基因組的不穩定性和突變是腫瘤的重要特征，癌基因和抑癌基因的異常是 ESCA 發生發展的重要原因[11]。目前，對 ESCA 致癌基因的研究顯著推動了其治療的進步。首先，是針對驅動基因改變的分子靶向治療。例如，EAC 傾向于過表達血管內皮生長因子（vascular endothelial growth factor，VEGF）形成新生血管發揮促腫瘤進展作用。針對 VEGF 信號通路的靶向治療已經成為晚期 EAC 的二線治療方案[12]。此外，EAC 可出現人類表皮生長因子受體-2 增殖相關通路的活化突變。無論是作為單一療法還是與 Abraxane 聯用，ramucirumab 已被美國食品藥品監督管理局（Food and Drug Administration，FDA）批準用于 EAC 的二線治療[13]。其次，針對免疫檢查點基因的免疫治療也是 ESCA 治療的方向和熱點。例如，針對程序性死亡分子 1（programmed death 1，PD-1）的單克隆抗體 pembrolizumab 已被 FDA 批準為聯合陽性評分≥10%的復發性、局部晚期或轉移性程序性死亡分子 1 配體（PD-1 ligand，PD-L1）陽性的 ESCC 的二線選擇。隨著組學技術的飛速發展和生物信息學技術的廣泛應用，建立基于基因標簽的臨床預測模型已成為構建新的腫瘤預后評估系統的趨勢和潮流。
本研究從癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）中獲取了 ESCC 和 EAC 的 RNA 測序數據和相應的臨床數據，建立了基于關鍵基因的 ESCC 和 EAC 的風險預測模型。隨后，我們對模型的準確性進行了評估。最后，我們分析了模型的臨床意義及其與腫瘤免疫微環境的關系。這些工作將為準確評估 ESCC 和 EAC 預后、優化治療手段提供重要依據。
1 資料與方法
1.1 數據資料
從 TCGA 數據庫（
1.2 統計學分析
利用“edgeR”包進行差異分析來獲取 ESCA 的差異基因［adjusted P value<0.05，|log2（fold change）| > 1］。通過單變量 Cox 回歸分析篩選出與總生存期（overall survival，OS）相關的差異基因。利用“glmnet”包進行最小絕對收縮和選擇算子（least absolute shrinkage and selection operator，LASSO）回歸分析，進一步明確 OS 相關核心基因。通過“survival”、“survminer”包進行多變量 Cox 回歸分析，建立 ESCC 和 EAC 的風險評分預后模型。以風險評分中位數為截斷值，將 ESCC 和 EAC 患者分別分成高、低風險組。利用“pheatmap”包繪制參與風險評分基因的熱圖。繪制 Kaplan-Meier 曲線并進行 log-rank 檢驗，探討風險評分與患者生存狀態的關系。通過“survival ROC”包進行受試者工作特征（receiver operating characteristic，ROC）曲線分析，評估基于風險評分的生存預測的靈敏度和特異度。ROC 曲線下面積代表了預后準確性。對臨床特征進行二分類，基于該分類的風險評分差異分析通過 Wilcoxon 檢驗實現。利用“rms”和“survival”包，構建了含臨床分期、T 分期、N 分期、風險評分參數的列線圖。通過校準曲線評估實際生存率和預測生存率的一致性。計算 C-index，評估預后預測模型的準確性。利用“estimate”包估計每個樣本在腫瘤微環境中的免疫成分和基質成分含量。基因集富集分析（gene set enrichment analysis，GSEA）闡述高風險組和低風險組之間差異表達的基因所富集的生物學過程。從 MSigDB 數據庫下載 KEGG 作為 GSEA 分析的參考基因集，只有 NOM P<0.05 的富集通路才被認為有意義。利用 CIBERSORT 估計腫瘤樣本的腫瘤浸潤免疫細胞豐度，僅 P<0.05 的腫瘤樣品被選擇出來。高、低風險組之間免疫細胞差異和免疫檢查點基因表達的差異分析均采用 Wilcoxon 檢驗。P≤0.05 為差異有統計學意義。所有統計學分析均利用 R 語言完成。
2 結果
2.1 ESCC 和 EAC 生存相關核心基因的識別
腫瘤差異基因在腫瘤發生發展中發揮重要作用，對腫瘤差異基因的識別有助于加深我們對腫瘤的認識。我們從 TCGA 數據庫獲取了 ESCA 的 RNA 測序數據，利用“edgeR”包得到了 ESCA 的差異基因（圖 1a）。這些差異基因包括 1 482 個上調基因和 1 765 個下調基因。隨后，我們在 ESCC 和 EAC 中對 ESCA 差異基因的表達進行了單變量 Cox 回歸分析。在 ESCC 中，180 個基因和 OS 相關（P<0.05）；在 EAC 中，371 個基因和 OS 相關（P<0.05，圖 1b）。考慮到我們得到的相關基因比較多，為進一步篩選關鍵基因，我們進行了 LASSO 回歸分析并計算了回歸系數（圖 1c～d）。在 ESCC 中，我們找到 3 個與 OS 顯著相關的基因（NUDCD1、NUP155、HLA-B），在 EAC 中，我們找到 6 個與 OS 顯著相關的基因（CIRBP、SLC26A9、MSH2、HIST1H3F、PMAIP1、HIST1H3J）。
圖1
ESCA 差異基因的獲取及 OS 相關基因識別
a：TCGA 數據庫中 ESCA 差異基因；b：單變量 Cox 回歸分析確定 ESCC 和 EAC 的 OS 相關基因；c：通過 LASSO 回歸分析進一步確定與 ESCC 的 OS 相關的 3 個關鍵基因；d：通過 LASSO 回歸分析進一步確定與 EAC 的 OS 相關的 6 個關鍵基因；TCGA：癌癥基因組圖譜；ESCA：食管癌；ESCC：食管鱗狀細胞癌；EAC：食管腺癌；OS：總生存期；LASSO：最小絕對收縮和選擇算子；FC：差異倍數；DEGs：差異表達基因
2.2 ESCC 和 EAC 風險預測模型的構建和評估
我們利用前期獲得的生存相關核心基因表達值進行了多變量 Cox 回歸分析（逐步回歸，表 1～2），并依托“survival”包中的“predict（）”函數算法得到了每個樣本的風險評分值。其中，ESCC 得到了由 3 個差異基因構成的風險評分預測模型，EAC 得到了由 5 個差異基因構成的風險評分預測模型。這些風險評分模型主要用于預測 ESCC 和 EAC 患者的預后。
表1
食管鱗狀細胞癌風險預測模型所涉及的基因
表2
食管腺癌風險預測模型所涉及的基因
以風險評分的中位數為截斷值，將 ESCC 和 EAC 患者分為高、低風險組。無論是 ESCC 患者還是 EAC 患者，高風險組患者存活人數更少（圖 2a～b）。
圖2
風險評分預后模型的構建和驗證
a：風險評分、患者生存時間的分布，ESCC 患者生存狀態和 ESCC 患者預后相關基因表達譜的熱圖；b：風險評分、患者生存時間的分布，EAC 患者生存狀態和 EAC 患者預后相關基因表達譜的熱圖；c：以風險評分為分層因素的 ESCC 患者的 Kaplan-Meier 生存分析；d：ROC 曲線分析顯示風險評分在預測 ESCC 患者生存時間中的準確性；e：以風險評分為分層因素的 EAC 患者的 Kaplan-Meier 生存分析；f：ROC 曲線分析顯示風險評分在預測 EAC 患者生存時間中的準確性；TCGA：癌癥基因組圖譜；ESCC：食管鱗狀細胞癌；EAC：食管腺癌；ROC：受試者工作特征；AUC：曲線下面積
風險評分在預測 ESCC 和 EAC 患者預后上所體現的價值需要進一步評估。單變量 Cox 回歸分析表明，ESCC 和 EAC 的風險評分均與患者 OS 顯著相關［HR=1.221，95%CI（1.116，1.336），P<0.001；HR=1.055，95%CI（1.023，1.087），P<0.001］。我們納入風險評分、TNM 分期和臨床分期進行了多變量 Cox 回歸分析。結果顯示，ESCC 和 EAC 的風險評分與 OS 顯著相關［HR=1.318，95%CI（1.173，1.482），P<0.001；HR=1.048，95%CI（1.016，1.082），P=0.003］。這表明，風險評分可作為 ESCC 和 EAC 患者的獨立預后因子。隨后，我們采用 Kaplan-Meier 法和 log-rank 檢驗進行生存分析，評估風險評分水平與患者 OS 的關系。結果顯示，高風險的 ESCC 和 EAC 患者 OS 顯著短于低風險組（P=0.009，P<0.001）。為評估風險預測模型的預測效果，我們進行了時間依賴的 ROC 曲線分析。ESCC 和 EAC 的風險評分預測1年生存率的曲線下面積（area under the curve，AUC）分別為 0.765 和 0.810，這提示風險預測模型預測效果較好（圖 2c～f）。
2.3 ESCC 和 EAC 風險預測模型的臨床意義
為進一步闡明風險評分與臨床特征的關系，我們將臨床特征進行二分類，通過 Wilcoxon 檢驗比較了這些臨床特征的風險評分差異。在 ESCC 中，T 分期為 T1～2 的患者風險評分更高（P=0.039），而性別（P=0.606）、淋巴結轉移與否（P=0.457）、有無遠處轉移（P=0.224）、腫瘤分期的高低（P=0.480）與風險評分的關系皆不顯著。在 EAC 中，淋巴結轉移（P=0.013）、存在遠處轉移（P=0.031）、臨床分期為 Ⅲ～Ⅳ期（P=0.016）患者風險評分更高，而 T 分期（P=0.965）、性別（P=0.214）與風險評分無關。
為預測 ESCC 和 EAC 患者的 OS，我們利用“rms”和“survival”包構建了包括風險評分和臨床特征（臨床分期、T 分期、N 分期）的列線圖。ESCC 和 EAC 患者列線圖的 C-index 分別為 0.870（SE=0.050）和 0.807（SE=0.046），展現出了較高的辨識能力。隨后，利用校準曲線對實際生存期和預測生存期之間的一致性進行評估。校準曲線顯示出了較好的預測能力（圖 3）。
圖3
ESCC 和 EAC 風險評分的臨床意義
a：以風險評分、臨床分期、T 分期、N 分期參數構建列線圖預測 ESCC 患者的 1 年生存率；b：通過校準曲線評價 ESCC 患者的預測 1 年生存率和實際 1 年生存率的一致性；c：以風險評分、臨床分期、T 分期、N 分期參數構建列線圖預測 EAC 患者的 1 年生存率；d：通過校準曲線評價 EAC 患者的預測 1 年生存率和實際 1 年生存率的一致性；TCGA：癌癥基因組圖譜；ESCC：食管鱗狀細胞癌；EAC：食管腺癌
2.4 風險預測模型與腫瘤微環境的關系
腫瘤微環境的結構成分主要由基質細胞和募集的免疫細胞構成。我們利用“estimate”算法對 ESCC 和 EAC 樣本中免疫成分和基質成分的數量進行了估計。為揭示腫瘤樣本中基質成分比例和免疫成分比例與風險評分的關系，采用 Wilcoxon 檢驗進行差異分析。結果表明，ESCC 中免疫成分比例更高的樣本風險評分更高（P=0.001），而基質成分比例高、低組風險評分差異無統計學意義（P=0.133）。在 EAC 樣本中，基質成分（P=0.387）和免疫成分（P=0.051）比例高、低組的風險評分差異均無統計學意義。
隨后，我們分別對高、低風險組的 ESCC 和 EAC 樣本進行了 GSEA 分析。結果顯示，低風險組的 ESCC 顯著富集在抗原呈遞和自然殺傷細胞介導的細胞毒性生物學過程（圖 4a～b）。這意味著風險評分越高的腫瘤部分免疫功能受到抑制。然而，基于風險評分分組的 EAC 與免疫相關的生物學過程關系并不密切，這與我們前面的結果吻合。
為進一步揭示風險評分與腫瘤浸潤免疫細胞的關系，我們利用 CIBERSORT 算法分析了 ESCA 組織中腫瘤浸潤免疫亞群的比例，并在 ESCC 和 EAC 樣本中分別構建了 21 種、20 種免疫細胞的表達譜（圖 4c～d）。隨后，我們評估了高低風險組的免疫細胞表達是否存在差異。結果表明，ESCC 高、低風險組樣本中活化的肥大細胞差異有統計學意義（P=0.035），EAC 高、低風險組樣品中免疫細胞差異無統計學意義（P>0.05，圖 4e～f）。隨后，我們探討了基于高、低風險分組的常見免疫檢查點基因的表達特征。我們發現，ESCC 患者中，大部分免疫檢查點基因在高風險組中高表達（P<0.05，圖 4g）。然而，在 EAC 患者中高、低風險組的免疫檢查點基因表達差異無統計學意義（P>0.05，圖 4h）。
3 討論
ESCA 是一種致命性惡性腫瘤，給家庭和社會帶來了沉重的負擔。ESCA 預后判斷的準確性影響著治療方法的選擇和患者預后的改善。因此，建立一套準確可行的預后判斷體系具有重要意義。腫瘤細胞基因表達的異常和當前組學技術的發展為建立基于基因的預后判斷模型提供了理論和技術支撐[11, 14]。
本研究中，我們利用 ESCA 腫瘤樣本和正常樣本的 RNA 測序數據進行差異分析，得到 1 482 個上調基因和 1 765 個下調基因。隨后，我們進行單變量 Cox 回歸分析，在這些差異基因中找到了 180 個與 ESCC 的 OS 相關基因，找到了 371 個與 EAC 的 OS 相關基因。由于篩選得到基因較多，我們進行了 LASSO 回歸分析、多變量 Cox 回歸分析，找到 3 個與 ESCC 的 OS 顯著相關基因（NUDCD1、NUP155、HLA-B），6 個與 EAC 的 OS 顯著相關基因（CIRBP、SLC26A9、MSH2、HIST1H3F、PMAIP1、HIST1H3J），并依托“survival”包中“predict（）”函數得到 ESCC 和 EAC 的風險預測模型。基于該風險預測模型的分組準確地將患者 OS 的長短區分出來。針對 1 年生存率的預測，ESCC 和 EAC 風險預測模型的 AUC 分別為 0.765 和 0.810，顯示出了較好的預測效果。此外，該風險模型與臨床特征關系密切。ESCC 的風險預測模型與腫瘤 T 分期相關，EAC 的風險預測模型與淋巴結轉移、遠處轉移和臨床分期相關。隨后，我們構建了包含風險預測模型和相關臨床特征的列線圖。ESCC 的風險預測模型與腫瘤樣本的免疫成分相關，這與我們的 GSEA 分析結果相吻合。通過進一步分析 ESCC 風險預測模型和免疫細胞的關系，我們發現 ESCC 的風險預測模型和肥大細胞激活途徑相關。我們以風險預測模型作為分組依據，評估了免疫檢查點基因的表達情況。我們發現，大部分免疫檢查點基因的表達和 ESCC 風險預測模型顯著相關。這進一步證實了免疫逃逸在腫瘤進展中的重要作用。
事實上，該評分模型所涉及的基因具有重要的生物學意義，尤其在腫瘤的轉化和惡性進展中。如NUDCD1是一種在多種腫瘤中高表達或激活的被稱為腫瘤抗原的癌蛋白[15-16]。NUDCD1 參與腫瘤細胞的增殖、侵襲、凋亡和上皮間質轉化等生物學過程，發揮促癌作用[16-17]。NUP155在肝細胞癌中高表達，與 TP53 突變和細胞周期調節基因的過表達相關，其高表達患者往往預后較差[18-19]。一項關于 ESCC 的研究[20]表明，NUP155 在 ESCC 中高表達，但是對 NUP155 進行的生存分析結果卻顯示高表達者預后更好。NUP155 在 ESCC 中具體作用并未闡明，仍需要進一步研究。HLA-Ⅰ類分子（HLA-A、HLA-B、HLA-C）的高表達與 ESCA、結直腸癌和肝細胞癌的較好預后有關，這可能歸因于腫瘤對 CD8+T 細胞抗腫瘤免疫力比較敏感[21-23]。而在胃癌、肺癌中，HLA-Ⅰ類分子高表達預示著更差的預后，可能是因為腫瘤從自然殺傷細胞中逃逸，這一過程被稱為“自我識別缺失”[24-25]。Zhang 等[26]認為，HLA-Ⅰ類分子表達的下調可能是 ESCC 的獨立不良預后因素。當然，有的研究[27]也認為 HLA-Ⅰ類分子表達與 ESCC 預后無關。Jiang 等[28]認為，HLA-B 在 ESCC 中上調，并在 ESCC 中發現了1個新的融合基因（HLA-E 和 HLA-B）。而我們的研究顯示，HLA-B 在 ESCC 中高表達，與預后不良相關。這與以前的研究存在一定的區別，這可能與腫瘤異質性或者腫瘤標本免疫細胞浸潤程度有關。CIRBP 屬于富含甘氨酸的 RNA 結合蛋白家族，在應激條件時發揮關鍵作用[29]。一方面，CIRBP 可以發揮抑癌作用，如抑制腫瘤細胞增殖[30-31]。另一方面，CIRBP 可發揮促癌作用，如促進腫瘤細胞增殖、減少凋亡和衰老、參與上皮間質轉化、增強侵襲和遷移等過程[29,32-35]。我們的研究表明，CIRBP 在 ESCA 中低表達，其表達值越高預示著 EAC 預后越好。CIRBP 在 EAC 中的作用暫不明確，需要深入探究。SLC26A9 主要作為 Cl?/HCO3?交換器、Cl?通道和 Na+耦合的轉運蛋白而發揮作用[36-37]。基于 TCGA 數據的分析表明，SLC26A9 在肺鱗狀細胞癌中低表達，但在基底乳腺癌中高表達[38]。我們的研究表明，SLC26A9 在 ESCA 中低表達，扮演保護基因的角色。但是，SLC26A9 在 EAC 中的具體作用還需進一步研究。MSH2 是第1個被鑒定的錯配修復（mismatch repair，MMR）基因，其蛋白功能是識別錯配的堿基[39-40]。MSH2 表達的缺失在 Barrett 相關 EAC 中較少見，它是部分 Barrett 相關 EAC 的 MMR 缺陷原因之一[41]。Evans 等[42]的研究表明，78% 的 EAC 中存在 MSH2 表達。我們的研究認為 MSH2 在 ESCA 中上調，這與之前的研究并不矛盾。HIST1H3F屬于依賴復制的組蛋白基因家族。HIST1H3F 在侵襲性乳腺癌中高表達，其高表達與無復發生存期、無遠處轉移生存期較差有關[43]。同樣，我們的研究認為 HIST1H3F 在 ESCA 中高表達，與預后不良有關。PMAIP1 是 Bcl-2 蛋白家族的促凋亡成員，其表達可激活線粒體凋亡級聯反應并誘導氧化應激和促進鈣釋放[44]。然而，PMAIP1 常在黑色素瘤中高表達，環磷腺苷效應元件結合蛋白介導的 MEK/ERK 信號通路激活可促進 PMAIP1 的表達以促進自噬。PMAIP1 介導的自噬激活在營養缺乏條件下對黑色素瘤細胞有間接的抗凋亡作用[45]。PMAIP1 在 EAC 中的作用需要進一步探索。
當然，我們的研究也存在一些缺陷。首先，我們的研究基于轉錄組數據，轉錄水平并不完全等同于蛋白水平。其次，轉錄組數據主要包含蛋白編碼基因，而非編碼 RNA 在腫瘤中亦發揮重要作用[46]。最后，我們的研究缺乏詳細的機制研究。這些工作將在以后的研究中展開。
利益沖突：無。
作者貢獻：馮楊波負責醞釀和設計實驗，實施研究，采集、分析和解釋數據，起草文章；熊延路負責醞釀和設計實驗，研究經費支持，分析和解釋數據，對文章內容作批評性審閱；趙晉波負責對文章內容作批評性審閱，統計學分析，指導研究；雷杰負責對文章內容作批評性審閱；辛少偉、喬天運負責采集數據，對文章內容作批評性審閱；周永生負責分析和解釋數據，對文章內容作批評性審閱；張曉負責統計分析指導，對文章內容作批評性審閱；姜濤負責對文章內容作批評性審閱；韓勇負責醞釀和設計實驗，研究經費支持，指導研究，對文章內容作批評性審閱及定稿。
