結(jié)腸癌淋巴結(jié)轉(zhuǎn)移的風(fēng)險基因及列線圖預(yù)測模型的構(gòu)建

武杰，李嵐，張惠博，吳思怡，宋啟斌

0 引言

結(jié)腸癌是常見的消化道惡性腫瘤之一，全球范圍內(nèi)結(jié)直腸癌發(fā)病率和死亡率分別居惡性腫瘤的第3位和第2位，結(jié)腸癌的5年總生存率約60%[1-3]。淋巴結(jié)轉(zhuǎn)移作為腫瘤分期的重要依據(jù)，對結(jié)腸癌患者預(yù)后的影響已被多項研究報道[4-6]。了解淋巴結(jié)轉(zhuǎn)移狀態(tài)可以為評估預(yù)后、指導(dǎo)治療等提供信息。手術(shù)病理診斷是評估淋巴結(jié)轉(zhuǎn)移的主要手段。NCCN指南指出術(shù)中清掃淋巴結(jié)12個以上才能保證淋巴結(jié)狀態(tài)的判斷準(zhǔn)確[7]；對于不能手術(shù)或淋巴結(jié)清掃不充分的患者，淋巴結(jié)狀態(tài)判斷可能不準(zhǔn)確。影像學(xué)也是了解淋巴結(jié)轉(zhuǎn)移的主要手段，但由于判斷的主觀性等限制，其判斷有一定偏倚，且其診斷敏感度可能弱于特異性[8]。

隨著基因測序技術(shù)的發(fā)展，從基因及分子的角度了解疾病更有助于推動疾病的精準(zhǔn)治療。腫瘤組織的基因表達(dá)值高低由客觀數(shù)值界定，避免了主觀偏倚。本研究旨在從基因?qū)用鎸ふ医Y(jié)腸癌淋巴結(jié)轉(zhuǎn)移的相關(guān)基因，并構(gòu)建由基因組成的nomogram預(yù)測模型，幫助評估患者的淋巴結(jié)狀態(tài)。

1 資料與方法

1.1 基因表達(dá)數(shù)據(jù)下載及過濾

從癌癥基因組圖譜 (The Cancer Genome Atlas，TCGA，https://portal.gdc.cancer.gov/) 數(shù)據(jù)庫下載獲得結(jié)腸癌患者的基因測序數(shù)據(jù)。納入標(biāo)準(zhǔn)：（1）非細(xì)胞株或動物來源樣本；（2）經(jīng)病理確診為結(jié)腸癌；（3）臨床信息包含明確的病理N分期；（4）診斷為N0期的患者需淋巴結(jié)清掃數(shù)目≥12個。不同探針對應(yīng)相同基因名時，取平均值進(jìn)行后續(xù)分析，并剔除低表達(dá)基因（在50%及以上的樣本中FPKM表達(dá)量為0的基因），以保證納入分析的基因有足夠的表達(dá)量。

1.2 差異表達(dá)基因分析及通路富集分析

經(jīng)以上標(biāo)準(zhǔn)過濾后，TCGA基因集被用于差異表達(dá)基因（differentially expressed gene,DEG）分析，分析結(jié)果中FDR（false positive rate）＜0.05且在淋巴結(jié)陽性和陰性組間表達(dá)均值差異大于兩倍（log2（fold change）＞1）的基因被認(rèn)為有差異表達(dá)。使用GSEA方法[9]對淋巴陽性和陰性組進(jìn)行基因集通路富集分析，鑒定兩組在GO[10]功能集及KEGG[11]上顯著富集的通路，富集顯著的標(biāo)準(zhǔn)設(shè)定為FDR＜0.05。

1.3 淋巴結(jié)轉(zhuǎn)移的風(fēng)險基因鑒定和模型建立

LASSO回歸[12]對差異表達(dá)基因進(jìn)一步降維，以AUC（area under receiver operating characteristic curve）值為標(biāo)準(zhǔn)，篩選轉(zhuǎn)移風(fēng)險關(guān)鍵基因。LASSO回歸基于L1正則化方法對樣本數(shù)據(jù)進(jìn)行變量選擇，將原本很小的系數(shù)直接壓縮至0，從而將這部分系數(shù)所對應(yīng)的變量視為非顯著性變量，將不顯著的變量直接舍棄。其原理主要基于L1正則化，是避免模型過擬合，將高維數(shù)據(jù)降維的有效方法。得到關(guān)鍵基因后利用赤池信息準(zhǔn)則（Akaike information criterion,AIC）構(gòu)建最優(yōu)的logistics列線圖（nomogram）模型。

1.4 模型的評估和驗證

采用ROC曲線及校準(zhǔn)曲線對模型的表現(xiàn)力進(jìn)行評估，同時對模型進(jìn)行Hosmer-Lemeshow擬合優(yōu)度檢驗。從NCBI平臺下的GEO 數(shù)據(jù)庫下載外部驗證數(shù)據(jù)。AUC（C統(tǒng)計量）大于0.65為理想模型。

1.5 決策曲線分析

決策曲線分析法[13-14]（decision curve analysis,DCA）用于探究在不同陽性閾值下，該模型對臨床凈獲益率的影響。DCA橫坐標(biāo)為閾概率（threshold probability）。當(dāng)nomogram模型評估值達(dá)到某值時，患者i的淋巴結(jié)轉(zhuǎn)移概率記為Pi；當(dāng)Pi達(dá)到某閾值（記為Pt），認(rèn)定為陽性。DCA的縱坐標(biāo)為凈獲益率（net benefit,NB）,定義為［真陽性數(shù)-假陽性數(shù)×Pt/(1-Pt)］/樣本數(shù)，由NB對Pt作圖可得DCA曲線。

1.6 統(tǒng)計學(xué)方法

差異表達(dá)基因分析采用非配對Wilcox檢驗，差異表達(dá)基因分析及富集分析中的P值經(jīng)過BH（Benjamin and Hochberg）法校正得到FDR值。模型變量的共線性通過方差膨脹因子（variance inflation factor,VIF）評估，VIF＞4時認(rèn)為變量間存在多重共線性。所有統(tǒng)計分析均在R軟件（https://www.r-project.org/）中實現(xiàn)，除差異表達(dá)分析及富集分析外，雙側(cè)P＜0.05為差異有統(tǒng)計學(xué)意義。

2 結(jié)果

2.1 患者基本信息

經(jīng)篩選，共納入TCGA數(shù)據(jù)庫中405例結(jié)腸癌患者的樣本，其中未發(fā)生淋巴結(jié)轉(zhuǎn)移的患者224例，發(fā)生淋巴結(jié)轉(zhuǎn)移的患者181例，基本情況見表1。在淋巴結(jié)陽性組中，＜60歲、男性、原發(fā)部位為左半結(jié)腸的患者較陰性組所占比例更高，淋巴結(jié)陰性組中，低T分期（Tis/T1/T2），低M分期（M0）及結(jié)腸息肉病史的患者所占比例更高。同時，淋巴結(jié)陽性組的CEA水平在統(tǒng)計學(xué)上顯著高于淋巴結(jié)陰性組患者。兩組淋巴結(jié)清掃數(shù)目無顯著差別。外部驗證數(shù)據(jù)來源于法國的多中心臨床數(shù)據(jù)集（儲存于GEO數(shù)據(jù)庫中，數(shù)據(jù)集編號GSE39582），其中未發(fā)生淋巴結(jié)轉(zhuǎn)移的患者312例，發(fā)生淋巴結(jié)轉(zhuǎn)移的患者241例。

2.2 差異表達(dá)基因分析及通路富集分析結(jié)果

差異表達(dá)分析共鑒定差異基因55個，其中陽性組上調(diào)基因46個，下調(diào)基因9個，見圖1A～B。GSEA富集分析得到與淋巴結(jié)轉(zhuǎn)移顯著正相關(guān)的GO通路551條，包括介導(dǎo)細(xì)胞間信號傳遞的表面受體通路、細(xì)胞發(fā)育、細(xì)胞突起形成、細(xì)胞膜形態(tài)改變等；負(fù)相關(guān)的GO通路45條，包括腫瘤細(xì)胞免疫、黏膜免疫、細(xì)胞殺傷等；正相關(guān)的KEGG富集通路28條，包括黏著斑、PI3K-Akt信號通路、細(xì)胞外基質(zhì)（ECM）受體相互作用、PPAR信號通路等；負(fù)相關(guān)KEGG通路38條，包括IL-17信號通路、Toll樣受體信號通路、PD1/PD-L1表達(dá)調(diào)控通路及抗原呈遞等，見圖1C～D。

表1 TCGA入選患者的臨床特征基線表Table 1 Clinical characteristics of included patients in TCGA database

2.3 LASSO回歸及風(fēng)險預(yù)測基因的鑒定

將編碼基因進(jìn)一步進(jìn)行LASSO回歸降維，共得到11個能有效預(yù)測淋巴結(jié)轉(zhuǎn)移的基因：TH、GREB1L、VWA5B1、PNMA6A、TNNC1、KIR2DL4、DRD1、STUM、SFTA2、CDH4和VGLL1，見圖2。對TCGA數(shù)據(jù)集中上述11個基因的高低表達(dá)進(jìn)行單因素及多因素logistic回歸分析得到淋巴結(jié)轉(zhuǎn)移的OR值，見表2。校正前后的OR值顯示CDH4、TH、GREB1L、PNMA6A、DRD1、VWA5B1、TNNC1、STUM、SFTA2的高表達(dá)提示淋巴結(jié)轉(zhuǎn)移的風(fēng)險增大，KIR2DL4的高表達(dá)提示淋巴結(jié)轉(zhuǎn)移風(fēng)險減小。VGLL1高表達(dá)在單因素分析中提示淋巴結(jié)轉(zhuǎn)移的風(fēng)險增大，而在多因素分析中提示淋巴結(jié)轉(zhuǎn)移的風(fēng)險減小。

表2 11個風(fēng)險基因的高低表達(dá)對淋巴結(jié)轉(zhuǎn)移狀態(tài)的OR值Table 2 Odd ratio values of 11 risk genes expression (high vs.low) for lymph nodes metastasis

2.4 基因nomogram預(yù)測模型的建立及驗證

根據(jù)逐步回歸分析的結(jié)果，由年齡、病理T分期、TH、CDH4、PNMA6A、TNNC1、KIR2DL4、STUM、SFTA2構(gòu)成的模型具有最小AIC（440.4）值。進(jìn)一步建立由上述變量組成的nomogram模型，見圖3，根據(jù)某患者每個變量實際情況找到其對應(yīng)的刻度，向上投射到頂部的標(biāo)尺（Points）即可得出每個變量的分值，分值相加即為總分值（total points），根據(jù)總得分值向下投射即可得到該患者的淋巴結(jié)轉(zhuǎn)移風(fēng)險概率（risk probability）。

得到的模型進(jìn)行內(nèi)部驗證發(fā)現(xiàn)模型的AUC值（AUC=0.800）、校準(zhǔn)度及擬合優(yōu)度（HL檢驗P＞0.05）均較佳，見圖4A～B，提示模型能很好地預(yù)測淋巴結(jié)轉(zhuǎn)移風(fēng)險。我們同時比較了納入或不納入基因的模型表現(xiàn)力，發(fā)現(xiàn)納入基因后模型的AUC值得到了較大的提高。在TCGA和外部驗證數(shù)據(jù)集中，AUC值的提高具有統(tǒng)計學(xué)意義（P＜0.001）。進(jìn)一步驗證發(fā)現(xiàn)，模型在法國多中心臨床數(shù)據(jù)集上的AUC值、校準(zhǔn)度及擬合優(yōu)度均較佳，說明風(fēng)險基因模型的外部數(shù)據(jù)適用性良好，見圖4C～D。

2.5 決策曲線分析及臨床影響曲線

決策曲線分析顯示模型在TCGA數(shù)據(jù)集及外部驗證數(shù)據(jù)集中，風(fēng)險閾值為0.2～0.7時，基于含風(fēng)險基因的nomogram模型進(jìn)行干預(yù)決策帶來的臨床獲益明顯高于未考慮風(fēng)險基因的模型，見圖4E～F。

圖1 差異表達(dá)基因及GESA功能富集分析Figure 1 Differentially expressed genes and GESA enrichment analysis

圖2 LASSO回歸法篩選關(guān)鍵基因Figure 2 Key genes screened by LASSO regression

3 討論

近年來，nomogram模型在癌癥領(lǐng)域中的應(yīng)用逐步增加，其將影響患者發(fā)病、預(yù)后或復(fù)發(fā)的各種臨床病理或基因等多個因素納入預(yù)測模型并可視化，將風(fēng)險比量化為具體分值，通過簡單運算即可獲得預(yù)測疾病復(fù)發(fā)、轉(zhuǎn)移及預(yù)后等的風(fēng)險概率，為臨床醫(yī)生和研究者提供了方便有利的工具[15-17]。本研究模型驗證結(jié)果提示其在預(yù)測淋巴結(jié)轉(zhuǎn)移狀態(tài)時有較好的一致性和區(qū)分度。GEO數(shù)據(jù)集上的驗證結(jié)果進(jìn)一步說明nomogram有良好的外部適用性。此外，決策曲線分析法表明使用nomogram評估可以在一定風(fēng)險閾值下帶來更高的臨床獲益。

圖3 Nomogram預(yù)測淋巴結(jié)轉(zhuǎn)移模型Figure 3 Nomogram model for predicting lymph node metastasis

淋巴結(jié)轉(zhuǎn)移對結(jié)腸癌預(yù)后的重要影響已被多項研究報道，伴有淋巴結(jié)轉(zhuǎn)移的患者較淋巴結(jié)陰性的患者生存更差，復(fù)發(fā)率也更高[18]。許多臨床病理因素可能對淋巴結(jié)轉(zhuǎn)移狀態(tài)有重要影響。Yamaok等[19]發(fā)現(xiàn)，隨著結(jié)腸癌浸潤深度的增加，淋巴結(jié)轉(zhuǎn)移率也逐漸升高。易小龍等[20]也發(fā)現(xiàn)腫瘤浸潤深度與分化程度是影響結(jié)腸癌淋巴結(jié)轉(zhuǎn)移的獨立因素。值得注意的是，在我們的研究中，＜60歲組患者的淋巴結(jié)轉(zhuǎn)移率高于≥60歲組，這可能與兩組人群具有不同的臨床特征有關(guān)[21]。影像學(xué)判斷結(jié)腸癌淋巴結(jié)轉(zhuǎn)移與否具有一定的局限性，通過聯(lián)合其他指標(biāo)可進(jìn)一步提高術(shù)前診斷淋巴結(jié)轉(zhuǎn)移的準(zhǔn)確性。因此，本研究的nomogram模型可進(jìn)一步協(xié)助臨床醫(yī)生判斷淋巴結(jié)轉(zhuǎn)移情況。

由于基因間可能存在很強(qiáng)的共線性，且考慮到結(jié)局事件數(shù)（淋巴結(jié)轉(zhuǎn)移數(shù)）遠(yuǎn)大于自變量數(shù)量，本研究先采用LASSO回歸進(jìn)行降維，尋找差異基因中與淋巴結(jié)轉(zhuǎn)移最為相關(guān)的基因。當(dāng)一組基因高度相關(guān)時，LASSO回歸會選出其中一個基因并且將其他因子收縮為零。LASSO回歸最終鑒定了11個與淋巴結(jié)轉(zhuǎn)移狀態(tài)高度相關(guān)的蛋白編碼基因。單、多因素分析提示僅KIR2DL4基因的高表達(dá)可能與抑制淋巴結(jié)轉(zhuǎn)移有關(guān)。KIR2DL4基因編碼殺傷細(xì)胞免疫球蛋白受體（killer-cell immunoglobulin-like receptor,KIR），在自然殺傷（natural killing,NK）細(xì)胞表面和少數(shù)T淋巴細(xì)胞表面表達(dá)，可識別并結(jié)合HLA-G分子，從而發(fā)揮對NK細(xì)胞的免疫抑制或激活的雙重作用[22]。值得注意的是，CDH4基因編碼的R-Cadherin蛋白（鈣黏蛋白家族成員之一），其功能在不同腫瘤中可能存在異質(zhì)性。據(jù)目前的研究報道，其可抑制鼻咽癌細(xì)胞的增殖和侵襲[23]，但在骨肉瘤和膠質(zhì)母細(xì)胞瘤中表現(xiàn)為促進(jìn)增殖和侵襲轉(zhuǎn)移的作用[24-25]。本研究中，CDH4的表達(dá)與結(jié)腸癌的淋巴結(jié)轉(zhuǎn)移呈正相關(guān)。

圖4 Nomogram模型的內(nèi)外部驗證Figure 4 I Internal and external validation of nomogram model

與以往研究相比，本研究具有一定的優(yōu)勢。首先，以往的研究僅僅關(guān)注了臨床病理因素[21,26]，這使淋巴結(jié)轉(zhuǎn)移的預(yù)測效果有限。而本研究中的模型加入基因后，模型的診斷預(yù)測能力顯著提高；其次，如前所述，影像學(xué)可以用于術(shù)前淋巴結(jié)轉(zhuǎn)移與否的判斷，但其主觀性較強(qiáng)，而本研究中的模型基于基因表達(dá)量進(jìn)行風(fēng)險評分，可以有效避免主觀判斷的偏倚；最后，本研究中的模型使用更加方便，nomogram模型將預(yù)測的風(fēng)險值量化，臨床工作者僅需通過簡單的計算即可獲得確定的淋巴結(jié)轉(zhuǎn)移概率值。但是，本研究也具有一定局限性，由于篇幅有限，對于11個風(fēng)險基因的具體生物學(xué)功能有待進(jìn)一步探究，且模型的適用性需在更大數(shù)據(jù)范圍內(nèi)進(jìn)行驗證。

綜上所述，本研究鑒定出了11個與淋巴結(jié)轉(zhuǎn)移相關(guān)的風(fēng)險基因，并構(gòu)建了由基因表達(dá)和臨床因素組成的nomogram預(yù)測模型。經(jīng)內(nèi)部驗證及外部驗證發(fā)現(xiàn)，該模型有良好的區(qū)分度和校準(zhǔn)度，并且可以帶來潛在的臨床獲益。