食管鱗癌淋巴結轉移關鍵基因鑒定及預后分析

李孟祥，程維剛，陳 攀，馮笑山，高社干，齊義軍

(河南科技大學 a.信息工程學院，b.臨床醫學院第一附屬醫院，中國 洛陽 471023)

食管癌是全球最常見的惡性腫瘤之一，發病率和死亡率分別位居所有惡性腫瘤的第七位和第六位[1]。2015年，中國食管癌的新發病例和死亡病例分別為24.6萬例和18.8萬例[2]，占全球新發病例和死亡病例的43%和36.9%。從組織學上看，中國食管癌90%以上為食管鱗狀細胞癌(Esophageal squamous cell carcinoma,ESCC)[3]。ESCC具有高度的侵襲和轉移能力，導致大部分ESCC患者初次確診時已發生腫瘤轉移，而腫瘤轉移是ESCC患者高死亡率、預后極差的主要原因之一。

食管黏膜及黏膜下層具有豐富的淋巴管網，并具有獨特的淋巴引流方式，極易發生食管內廣泛或跳躍性轉移及頸、胸、腹三野淋巴結轉移。淋巴結轉移與腫瘤浸潤深度和腫瘤分化程度密切相關，術中淋巴結清掃數目、淋巴結轉移陽性數目、淋巴結轉移陽性率等是判定ESCC預后的重要獨立危險因素[4-6]。此外，淋巴結轉移也是ESCC TNM分期中重要的病理特征之一，直接影響ESCC患者術后治療方案的選擇[7]。因此，通過基因差異表達譜分析，確定與淋巴結轉移高度相關的分子變異，對于ESCC個體化精準治療極其重要。

極端梯度提升(Extreme gradient boosting,XGBoost)是基于Boosting集成的一種機器學習算法，適用于大規模數據的分布式并行運算。該算法具有高維數據處理、分析缺失值、運算效率高及可遷移性強等優點，近年來被廣泛應用于數據挖掘。本研究利用XGBoost算法分析基因表達數據庫(Gene expression omnibus,GEO)中ESCC mRNA轉錄組數據，鑒定ESCC淋巴結轉移關鍵mRNA分子譜，構建淋巴結轉移分類模型。

1 材料與方法

1.1 ESCC mRNA表達譜數據

從GEO下載數據集GSE53624和GSE53622，這兩個數據集分別包括119例和60例ESCC和配對癌旁組織基因表達譜數據，基因芯片平臺為GPL18109(Agilent-038314 CBC Homo sapiens lncRNA + mRNA microarray V2.0)。從https://www.agilent.com/下載平臺GPL18109的探針組序列，利用GENCODE和SeqMap進行序列比對，重新注釋數據，提取mRNA表達譜數據。

1.2 ESCC淋巴結轉移相關mRNA鑒定

根據文獻報道方法[8]，為減少不同樣本之間異質性的影響，將ESCC與配對癌旁組織表達值的差值作為ESCC mRNA表達譜數據。根據淋巴結轉移與否，將ESCC樣本分為淋巴結轉移陰性組(N0組)和淋巴結轉移陽性組(N1組，包括臨床分期中N1，N2和N3)，篩選淋巴結轉移相關mRNA分子。將GSE53624中119例ESCC樣本隨機分為訓練集(60例)和測試集(59例)，GSE53622中60例ESCC樣本作為獨立驗證組。表1顯示了3個數據集中ESCC樣本的人口統計學和臨床病理特征，包括年齡、性別、腫瘤部位、腫瘤分級、T stage及TNM stage等。在60例訓練集中，以student’s T檢驗P<0.05和淋巴結轉移陽性組與陰性組間表達值之差大于0.5為篩選標準，鑒定ESCC淋巴結轉移相關的差異表達mRNA分子。

1.3 算法簡介

XGBoost算法以CART分類樹為基學習器，來源于Boosting方法。在迭代過程中，后一個模型對前一個模型的誤差進行校正，通過擬合殘差優化目標函數，提高預測分類準確率。XGBoost對損失函數進行二階泰勒展開，并在損失函數中加入正則項以控制目標函數的下降和模型復雜度，防止模型過擬合[9]。邏輯回歸(Logistic regression,LR)是一種廣義的線性模型，是在線性回歸的基礎上外加一層Sigmoid函數映射。支持向量機(Support vector machine,SVM)算法通過非線性映射將輸入空間映射到一個高維空間，進而在高維空間中構造最優分類超平面，利用支持向量最大化幾何間隔，降低分類誤差。本研究分別使用R 語言中bestglm，e1701和xgboost等程序包實現LR，SVM和XGBoost模型構建。

1.4 模型預測效能評價

本文應用R 3.63進行統計學分析。使用受試者操作特征(Receiver operating characteristic,ROC)曲線比較每個預測模型的敏感性和特異性，并以ROC曲線下面積(Area under ROC curve,AUC)值評價作為預測結果，AUC值反映分類模型預測的精準度。所有ROC曲線應用pROC函數包進行計算。K-S(Kolmogorov-Smirnov)值是另一種分類模型效能的評價指標，K-S曲線將選定的閾值作為橫軸，分類模型的真陽性率(True positive rate,TPR)和假陽性率(False positive rate,FPR)均繪制到縱軸，K-S值為所有TPR和FPR差值的絕對值中的最大值，K-S值的大小與分類模型區分特征的準確性呈正相關。

1.5 生物學功能及通路富集分析

對ESCC淋巴結轉移關鍵基因進行Gene Ontology(GO)富集分析，包括細胞成分(Cellular component,CC)、生物學過程(Biological process,BP)和分子功能(Molecular function,MF)，根據超幾何分布檢驗的錯誤發現率(False Discovery Rate,FDR)，確定關鍵基因的細胞定位、分子功能及參與的生物學過程。使用ClusterProfiler程序包進行富集分析及相關繪圖[10]。

1.6 生存分析

采用Kaplan-Meier法繪制生存曲線，對數秩和檢驗(Log-Rank test)進行生存期差異顯著性檢驗。連續變量用survminer包中函數surv_cutpoint確定最佳截斷值，將連續變量簡化為二分類變量。再用單因素和多因素Cox比例風險回歸模型確定預后的影響因素，似然比檢驗(Likelihood ratio test)確定模型顯著性，并根據多變量分析結果繪制森林圖。應用survival函數包進行生存分析。

2 結果

2.1 ESCC淋巴結轉移相關差異表達基因

分析GSE53624訓練集中60例ESCC樣本mRNA表達譜數據，根據t檢驗P<0.05且兩組之間mRNA表達值均數之差>0.5為篩選條件，鑒定了ESCC淋巴結轉移相關的159個差異表達基因，包括淋巴結轉移陽性組中31個高表達和128個低表達的mRNA分子(圖1a)。

2.2 ESCC淋巴結轉移關鍵分子

以上述的159個ESCC淋巴結轉移相關mRNA分子作為XGBoost模型的初始特征集合，網格搜索和5倍交叉驗證方法確定XGBoost模型超參數，包括最大迭代次數(nrounds=200)、學習率(eta=0.1)、單棵樹最大深度(max_depth=4)、最小減損函數下降值(gamma=0.2)、隨機采樣特征比率(colsample_bytree=0.3)以及葉子節點最小權重(min_child_weight=0.7)等，用GSE53624訓練集60例ESCC樣本建立ESCC淋巴結轉移預測模型。根據Gain值評價159個mRNA分子重要性，圖1b顯示了159個差異表達mRNA的重要性分布，其中18個mRNA分子的重要性分值大于0.02，其余mRNA分子的重要性分值小于0.02。因此，本研究將重要性分值>0.02的18個mRNA分子作為ESCC淋巴結轉移的關鍵mRNA分子，將其納入分類模型。18個mRNA分子及其Gain值分別為：MASP1(0.068 5)，ANOS1(0.063 2)，CENPP(0.058 5)，ABCG2(0.043 8)，GALNT12(0.035 1)，IP6K3(0.034 5)，SLC16A5(0.032 6)，MMP27(0.030 7)，C6orf15(0.030 5)，KRT6C(0.028 6)，CXCL10(0.028 1)，RIMS2(0.027 3)，RPTN(0.026)，LIMA1(0.025 1)，KRT6B(0.024 9)，TNC(0.022 6)，LCE3D(0.021 8)和APLF(0.021 5)。

2.3 構建ESCC淋巴結轉移分類模型

為構建ESCC淋巴結轉移分類器，從重要性分值最大的兩個mRNA分子開始依次增加mRNA分子，構成逐漸增大的mRNA分子集合，用XGBoost，LR和SVM算法計算2～20個mRNA組合在測試集中預測ESCC淋巴結轉移的AUC值。在XGBoost算法中，重要性分值最高的兩個mRNA分子模型在測試集中診斷ESCC淋巴結轉移的AUC值為0.618 3，隨著特征性mRNA分子增加，分類模型的診斷效能逐漸升高，18個mRNA分子模型的AUC值最大(0.793 7，圖1c)；而LR與SVM算法在特征分子增加過程中，最大的AUC值分別為0.713和0.714，預測結果不穩定(圖1c)。

為比較本研究建立的18個mRNA分子預測ESCC淋巴結轉移的效能(XGB-18 mRNA)，應用同樣的18個mRNA分子，建立了LR-18 mRNA和SVM-18 mRNA的分類模型。LR-18 mRNA和SVM-18 mRNA在同一的訓練集上進行超參數及核函數選擇。XGB-18 mRNA，LR-18 mRNA和SVM-18 mRNA 3種模型在測試集中診斷ESCC淋巴結轉移的AUC值分別為0.793 7，0.676 0和0.695 8，XGB-18 mRNA模型預測ESCC淋巴結轉移的效能顯著高于LR模型(Z=2.35,P=0.018)和SVM模型(Z=1.60,P=0.11)，其ROC曲線見圖1d。在外部驗證集(GSE53622)中，XGB-18 mRNA，LR-18 mRNA和SVM-18 mRNA 3個模型的ROC曲線如圖1e所示，AUC值分別為0.711，0.669和0.673，可見XGB-18 mRNA模型的AUC值最高。此外，利用模型區分度評價指標K-S值對XGB-18 mRNA，LR-18 mRNA和SVM-18 mRNA 3個模型進行評估(圖1f-h)，其K-S值分別為0.468，0.242和0.300，表明XGB-18 mRNA模型預測ESCC淋巴結轉移的準確性最高。

注：紅色和綠色分別代表在淋巴結轉移陽性組中表達量高于和低于淋巴結轉移陰性組的mRNA。圖1 特征mRNA選擇、模型建立與模型效果評價(a)訓練集中差異表達mRNA火山圖；(b)XGBoost模型中159個mRNA的重要性評分；(c)不同mRNA特征子集的預測性能；(d)3種預測模型在測試集上的ROC曲線比較；(e)3種模型在外部驗證集中的ROC曲線比較；(f-h)3個模型的K-S(Kolmogorov-Smirnov)值Fig. 1 Feature selection,model construction and prediction efficiency (a)Volcano plot for the differential expression of mRNAs in the training set，(b)the importance scores of 159 mRNAs assessed by the XGBoost algorithm,(c)prediction performances of different mRNA feature subsets,(d and e)receiver operating characteristic curves of three prediction models in the testing set and the external validation set,and (f-h)Kolmogorov Smirnov values of the three models

2.4 生存分析

為明確XGB-18 mRNA分類模型的臨床意義，本實驗分析了XGB-18 mRNA分類模型在ESCC生存預后中的作用。在59例測試集和60例外部驗證集中，以XGB-18 mRNA預測值的最佳界值將樣本分為高風險組和低風險組。XGB-18 mRNA高低風險組ESCC患者生存曲線表明，XGB-18 mRNA高風險組ESCC患者的總體生存時間低于低風險組，GSE53624測試集中XGB-18 mRNA高、低風險組ESCC患者的中位生存時間分別為12.6和56.2月，生存分析顯示該模型預測值為預后危險因素(HR 3.91,95%CI 1.95～7.84;P<0.000 1，圖2a)。進一步做亞組分析，淋巴結轉移陰性患者中XGB-18 mRNA高風險組ESCC患者生存時間也明顯低于低風險組(HR 6.75,95%CI 1.59～28.76;P=0.003 3，圖2b)，淋巴結轉移陽性患者中，XGB-18 mRNA高風險組ESCC患者生存時間明顯較低風險組短(HR 3.11,95%CI 1.24～7.82;P=0.012，圖2c)。在外部驗證集中，高、低風險組中位生存時間分別為10.2和39.8月，生存分析顯示模型預測值為預后風險因素(HR 2.27,95%CI 0.99～5.25;P=0.048 9，圖2d)；亞組分析結果與測試集中類似。

圖2 測試集和驗證集ESCC患者的生存分析(a)測試集59例ESCC患者中，高和低XGB-18 mRNA score組的生存曲線；(b-c)測試集中分別以淋巴結轉移陰性、陽性分組做亞組分析，高和低XGB-18 mRNA score組的生存曲線；(d)外部驗證集60例ESCC中，高和低XGB-18 mRNA score組的生存曲線；(e-f)驗證集中分別以淋巴結轉移陽性、陰性分組做亞組分析，高和低XGB-18 mRNA score組的生存曲線Fig. 2 Survival analysis of ESCC in the testing set and validation set (Survival curves of 59 ESCC patients with high-and low-XGB-18 mRNA scores in the testing set (a),and the subsets of N0 (b)and N1 (c);survival curves of 60 ESCC patients with high-and low-XGB-18 mRNA scores in the validation set (d),and the subsets of N0 (e)and N1 (f))

2.5 單因素和多因素Cox回歸模型分析

將測試集和驗證集中ESCC患者的人口學、臨床病理特征及XGB-18 mRNA模型預測值score進行單因素Cox回歸模型分析，發現N stage、TNM stage和XGB-18 mRNA預測值是ESCC預后影響因子(P<0.05，表2)。多因素Cox回歸模型分析結果表明，XGB-18 mRNA預測值score(HR 3.3,95%CI 1.48～7.6;P=0.004，表3)是影響ESCC預后的獨立因子。在GSE53622驗證集中，多因素Cox回歸模型分析結果與測試數據集類似(HR 4.0,95%CI 1.52～10.7;P=0.005，表3)。

表2 測試集和驗證集ESCC患者的人口學特征、臨床病理特征和XGB-18 mRNA預測值score的單因素Cox回歸模型分析Tab. 2 Univariate Cox regression analyses of the demographics characteristics,clinicopathological characteristics and XGBoost-18 mRNA scores of ESCC from the test set and the validation set

表3 測試集和驗證集上ESCC患者的多因素Cox回歸模型分析Tab. 3 Multivariate Cox regression analyses of ESCC from the test set and the validation set

2.6 GO富集分析

對XGB-18 mRNA模型納入的18個mRNA 進行GO富集分析，在生物學過程中，上皮角化過程顯著富集，參與該生物學過程的有4個mRNA分子，分別是KRT6B，KRT6C，LCE3D和RPTN等(圖3a)。這4個mRNA分子在ESCC中均為低表達，并且在淋巴結轉移陽性ESCC中的表達進一步降低，提示淋巴結轉移陽性的ESCC癌細胞角化過程受阻，使癌細胞不同程度地失去上皮特征性分子表達，可能發生了上皮間質轉化(圖3d-g)。細胞定位富集分析結果表明，18個mRNA分子主要位于細胞內不溶性膜、角質微絲、突觸周圍和突觸相關細胞骨架(圖3b)。分子功能富集分析結果表明，18個mRNA分子的功能主要包括己烷基磷酸肌醇6，1，3，5等激酶活性、CXCR3受體結合、異生型跨膜ATP酶活性和cAMP依賴的蛋白激酶調節活性(圖3c)。

3 討論

由于食管黏膜及黏膜下層具有豐富的淋巴管網，因此淋巴結轉移是食管癌轉移的主要方式，浸潤至黏膜下層的食管癌細胞發生淋巴結轉移陽性率可達20%～30%，而局限于黏膜層的食管癌較少發生淋巴結轉移。淋巴結轉移是食管癌進展期的主要特征，與食管癌患者的術后生存預后密切相關，多學科綜合治療是進展期食管癌最佳治療方案，而預后預測及分型是食管癌患者進行個體化治療的關鍵，故而精準診斷ESCC淋巴結轉移是ESCC臨床治療的關鍵[11,12]。

本研究基于ESCC mRNA表達譜數據，通過差異分析確定了與ESCC淋巴結轉移相關的159個mRNA分子，并將其作為輸入特征結合機器學習算法建立模型。根據這些mRNA分子特征在模型中的重要性，最終建立了由18個mRNA分子組成的XGB-18 mRNA ESCC淋巴結轉移診斷模型。18個mRNA分子的生物學功能富集分析表明(圖3a-c)，發生淋巴結轉移的ESCC細胞角質化過程受阻，KRT6B，KRT6C，LCE3D及RPTN等與細胞角質化相關的4個基因在淋巴結轉移陽性ESCC中表達均下調。正常食管鱗狀上皮細胞表達上皮細胞粘附分子(Ep-CAM)、鈣粘附蛋白(E-cadherin)、緊密連接蛋白(Occludin，ZO-1)和細胞骨架相關(Keratin，ezrin),維系上皮組織結構完整性，防止物理、化學、微生物、炎癥分子等有害因子的損傷作用。在食管上皮癌變過程中，癌變細胞失去這些結構性分子，使細胞間粘附力降低，獲得遷移和侵襲能力。由此可見，參與細胞角質化過程的基因可能作為淋巴結轉移陽性ESCC的診斷和臨床治療的靶分子。

圖3 18個ESCC淋巴結轉移關鍵分子的GO富集分析與4個角質化相關分子在ESCC和癌旁組織中的表達(a)生物過程(BP)富集分析結果；(b)細胞組成(CC)富集分析結果；(c)分子功能(MF)富集分析結果；(d-g)分子KRT6B，KRT6C，LCE3D和RPTN在淋巴結轉移陰性組(N0)、淋巴結轉移陽性組(N1)以及所有179例樣本中的表達量

由于病變部位、手術難度等因素影響，不同ESCC患者在手術過程中淋巴結清掃數目和區域存在較大差異。癌癥基因圖譜數據庫(The Cancer genome atlas,TCGA)中，92例ESCC中包括53例ESCC患者進行了手術淋巴結清掃，淋巴結清掃數目從1到83個不等。2015年鄭州大學附屬腫瘤醫院胸外科的一項研究表明[13]，2010—2014年治療的313例胸段ESCC患者中有122例(38.97%)發生了淋巴結轉移，313例患者共清掃淋巴結4 461枚(平均14.2枚/人)，其中癌轉移的淋巴結有294枚，淋巴結轉移率為6.59%。其中，喉返神經旁淋巴結轉移率最高(25.5%)，賁門胃左動脈旁次之(18.2%)。另一項南京醫科大學一附院腫瘤科分析了1 791例食管癌樣本淋巴結轉移情況[14]，其中1 693例ESCC(占94.5%)樣本中發生淋巴結轉移的患者有586例(34.61%)，術中共清掃淋巴結17 674枚(平均10.4枚/人)，淋巴結轉移陽性為1 409枚，淋巴結轉移率為7.97%。以上研究表明，ESCC術中淋巴結清掃數目和清掃區域具有較大的異質性，淋巴結清掃數目過少、陽性淋巴結未被檢測出等導致淋巴結轉移陰性患者可能存在假陰性。2017年關于食管癌根治術胸部淋巴結清掃的中國專家共識中指出，2016版NCCN建議的二野或三野淋巴結清掃數目須達到11～15枚[15]，才能提供較為準確的N分期和TNM分期判斷，正確指導術后治療。此外，本研究包含的病例資料中未提供淋巴結轉移區域信息，可能導致預后預測結果存在誤差。

目前，XGBoost算法在生物醫學領域應用廣泛，據此建立的生物學模型的診斷效果較為理想[16-18]。有研究表明[16]，經過數據歸一化、超參數選優后，用XGBoost算法以肝功能、血脂、腎功能、乙肝、血常規指標等為特征建立血糖值的回歸模型，以均方根誤差和平均絕對百分比誤差作為模型評價指標，結果表明基于XGBoost算法的模型具有精度高、運行快、穩定性強等優勢，較基于SVM的模型和隨機森林模型預測精確。另一研究[17]利用XGBoost算法挖掘與N2-3期淋巴結轉移相關的因素，通過交叉驗證和多次隨機化實驗得到最佳變量組合，該組合包括腫瘤大小、組織學類型、多灶性、淋巴管浸潤、ER陽性細胞百分比和前哨淋巴結陽性數目等，結果表明，基于XGBoost算法的預測模型優于邏輯回歸算法模型，達到0.80 (95%CI 0.65～0.92)。

總之，本文利用XGBoost算法鑒定了18個ESCC淋巴結轉移關鍵基因，并建立了XGB-18 mRNA分類模型，其預測ESCC淋巴結轉移的準確性高于LR模型和SVM模型，并且XGB-18 mRNA模型預測的風險值，是ESCC預后獨立危險因素，為ESCC患者臨床個體化治療策略及方案制定提供理論依據和指導。