包汝娟，陳慧芳，董 宇，葉幼瓊，蘇 冰

上海交通大學(xué)醫(yī)學(xué)院上海市免疫學(xué)研究所，上海200025

結(jié)直腸癌（colorectal cancer，CRC）是最常見的惡性腫瘤之一，發(fā)病率和死亡率均較高［1］。目前，CRC 患者的預(yù)后主要取決于診斷時(shí)的疾病進(jìn)程和腫瘤分期［2］。然而，由于缺乏適當(dāng)?shù)脑\斷方法，預(yù)后指標(biāo)并不能滿足實(shí)際的臨床需求［3］。因此，開發(fā)能夠有效預(yù)測(cè)患者預(yù)后的評(píng)價(jià)指標(biāo)十分必要。目前，以靶向程序性死亡受體-1（programmed death-1，PD-1）為代表的免疫檢查點(diǎn)抑制劑是CRC 免疫治療的主要研究方向，但僅有一小部分CRC 患者受益，且在其他癌癥中免疫治療也并非適合于所有患者［4-5］。因此，找到能夠預(yù)測(cè)癌癥患者免疫治療效果的關(guān)鍵特征基因，將有助于指導(dǎo)癌癥的臨床治療。近年來，隨著微陣列技術(shù)［6］和RNA 測(cè)序技術(shù)的飛速發(fā)展，由此產(chǎn)生的轉(zhuǎn)錄組數(shù)據(jù)成為了探索疾病機(jī)制、挖掘各類疾病預(yù)后標(biāo)志物的豐富來源［7-9］。因此，整合來源于不同技術(shù)平臺(tái)的數(shù)據(jù)集，將會(huì)為全轉(zhuǎn)錄組水平分析腫瘤相關(guān)異?；蛱峁└C合的數(shù)據(jù)基礎(chǔ)。

因此，本研究利用生物信息學(xué)技術(shù)整合公共數(shù)據(jù)庫的CRC 數(shù)據(jù)集，對(duì)CRC 預(yù)后相關(guān)的關(guān)鍵基因進(jìn)行篩選并構(gòu)建模型，評(píng)估模型的預(yù)測(cè)性能，探索模型在預(yù)測(cè)免疫治療效果方面的應(yīng)用，為CRC 患者的預(yù)后管理和個(gè)體化精準(zhǔn)免疫治療提供參考。

1 資料與方法

1.1 CRC相關(guān)數(shù)據(jù)的獲取及整理

RNA-seq 數(shù)據(jù)集：①從基因型-組織表達(dá)（genotypetissue expression，GTEx）數(shù)據(jù)庫（https：//commonfund.nih.gov/gtex）中下載308例正常結(jié)直腸組織樣本的基因表達(dá)數(shù)據(jù)及對(duì)應(yīng)的捐贈(zèng)者臨床信息（包含捐贈(zèng)者的性別、組織來源等），即為GTEx 數(shù)據(jù)集，將基因表達(dá)數(shù)據(jù)中每千個(gè)堿基的轉(zhuǎn)錄每百萬映射讀取的片段數(shù)（fragments per kilobase of exon model per million mapped fragments，F(xiàn)PKM）值進(jìn)行l(wèi)og2（x+0.001）轉(zhuǎn)換。②在癌癥基因組圖譜（the cancer genome atlas，TCGA） 數(shù)據(jù)庫（https：//portal.gdc.cancer.gov）中下載471例CRC組織樣本及41例正常結(jié)直腸組織樣本的基因表達(dá)數(shù)據(jù)以及對(duì)應(yīng)的患者臨床信息（包含患者的年齡、性別、生存狀態(tài)、生存時(shí)間等），即為TCGA 數(shù)據(jù)集，對(duì)基因表達(dá)數(shù)據(jù)的FPKM 值進(jìn)行l(wèi)og2（x+1）轉(zhuǎn)換。隨后，使用sva R包的ComBat函數(shù)將GTEx 數(shù)據(jù)集和TCGA 數(shù)據(jù)集整合為一個(gè)數(shù)據(jù)集，即為RNA-seq數(shù)據(jù)集。

微陣列數(shù)據(jù)集：在基因表達(dá)綜合（gene expression omnibus，GEO） 數(shù) 據(jù) 庫［10］（https：//www.ncbi.nlm.nih.gov/geo/query/acc.cgi）中下載10 例CRC 患者的基因表達(dá)數(shù)據(jù)及對(duì)應(yīng)的患者臨床信息（包含患者的年齡、性別、生存狀態(tài)、生存時(shí)間等），包括GSE8671［11］、GSE18105［12］、GSE20916［13］、GSE23878［14］、GSE37364［15］、GSE21510［16］、GSE33113［17］、GSE39582［18］、GSE17536［19］和GSE17537［19］。使 用oligo R 包（http：//www.bioconductor.org/packages/release/bioc/html/oligo.html）的Robust Multiarray Averaging函數(shù)對(duì)原始數(shù)據(jù)進(jìn)行處理。由于前8個(gè)數(shù)據(jù)集同時(shí)包含了CRC 組織和正常結(jié)直腸組織的樣本，故用于后續(xù)差異分析；后3 個(gè)數(shù)據(jù)集包含了模型研究部分所需的臨床信息，故而作為模型驗(yàn)證集進(jìn)行后續(xù)研究。

1.2 差異表達(dá)基因篩選及功能分析

使用R 語言limma 包［20］分別對(duì)每個(gè)微陣列數(shù)據(jù)集中的差異表達(dá)基因（differentially expressed genes，DEGs）進(jìn)行篩選，再通過R 語言RobustRankAggreg （RRA）包［21］獲得前8 個(gè)微陣列數(shù)據(jù)集中共有的DEGs。基因表達(dá)的差異用P 值和差異倍數(shù)（fold change，F(xiàn)C）的對(duì)數(shù)（log2FC） 表 示。將P＜0.05 且|log2FC|＞1 的 基 因 視 為DEGs。隨后，分別對(duì)RNA-seq 數(shù)據(jù)集、TCGA 數(shù)據(jù)集中的DEGs 進(jìn)行篩選，方法及篩選標(biāo)準(zhǔn)同上。使用R 語言clusterProfiler 包［22］分別對(duì)微陣列數(shù)據(jù)集、RNA-seq 數(shù)據(jù)集獲得的DEGs 進(jìn)行基因本體數(shù)據(jù)庫（Gene Ontology，GO）功能分析，結(jié)果以P＜0.01 為入選標(biāo)準(zhǔn)。采用R 語言GOplot包［23］呈現(xiàn)GO功能分析的結(jié)果。

1.3 預(yù)后風(fēng)險(xiǎn)評(píng)分模型的構(gòu)建和性能評(píng)估

1.3.1 用于構(gòu)建及評(píng)估模型的數(shù)據(jù)集說明 因建模構(gòu)建需要，結(jié)合數(shù)據(jù)集的基因表達(dá)數(shù)據(jù)及臨床信息中的生存狀態(tài)、生存時(shí)間，從TCGA數(shù)據(jù)集中選擇了438例符合上述要求的樣本開展研究，其中隨機(jī)抽取219例樣本作為訓(xùn)練集，排除訓(xùn)練集的剩余219例樣本作為內(nèi)部驗(yàn)證集；同時(shí)，從微陣列數(shù)據(jù)集中選擇GSE39582（556 例樣本）、GSE17536（177 例樣本）和GSE17537（55 例樣本）數(shù)據(jù)集作為外部驗(yàn)證集。

1.3.2 預(yù)后風(fēng)險(xiǎn)評(píng)分模型的構(gòu)建及分組標(biāo)準(zhǔn) 對(duì)“1.2”中已篩選獲得的DEGs 進(jìn)行二次篩選，獲得GSE39582、TCGA 數(shù)據(jù)集這2 個(gè)數(shù)據(jù)共有的DEGs。采用R 語言survival包［24］進(jìn)行單因素Cox回歸分析，篩選與不良預(yù)后相關(guān)的基因［風(fēng)險(xiǎn)比（hazard ration，HR）＞1，P＜0.05］，并將其作為構(gòu)建模型的候選基因?；谟?xùn)練集，通過R語言glmnet 包［25］將候選基因作為參數(shù)進(jìn)行套索（LASSO）回歸分析，確定最佳懲罰值，而后選擇其對(duì)應(yīng)的回歸系數(shù)不為0 的基因作為建?；颍傩卸嘁蛩谻ox 回歸分析，構(gòu)建預(yù)后風(fēng)險(xiǎn)評(píng)分模型。風(fēng)險(xiǎn)評(píng)分計(jì)算公式為：風(fēng)險(xiǎn)評(píng)分=基因1表達(dá)量×多因素Cox回歸系數(shù)1+…+基因N表達(dá)量×多因素Cox回歸系數(shù)N（其中N代表基因數(shù)）。計(jì)算每例患者的風(fēng)險(xiǎn)評(píng)分，通過R 語言survival 包［24］的surv_cutpoint函數(shù)將患者劃分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。

1.3.3 在訓(xùn)練集中評(píng)估預(yù)后風(fēng)險(xiǎn)評(píng)分模型的性能 采用上述計(jì)算公式對(duì)訓(xùn)練集中每例患者的風(fēng)險(xiǎn)評(píng)分進(jìn)行計(jì)算，并據(jù)此將其分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。通過R語言ggplot2包繪制2組患者的風(fēng)險(xiǎn)評(píng)分和生存狀態(tài)的分布圖，以及建?；虮磉_(dá)量圖譜。通過R語言survivalROC包［26］繪制時(shí)間依賴性受試者操作特征曲線（receiver operator characteristic curve，ROC 曲線），計(jì)算訓(xùn)練集中患者生存時(shí)間分別在1、2、3年的曲線下面積（area under the curve，AUC）；通過R語言survminer包（https：//cran.r-project.org/web/packages/survminer/index. html） 繪 制Kaplan-Meier（KM）生存曲線，計(jì)算2 組患者的生存率。根據(jù)AUC 值和組間的KM 生存率的差異，評(píng)估模型在訓(xùn)練集中的預(yù)測(cè)性能。同時(shí)，結(jié)合訓(xùn)練集中CRC 患者的年齡、性別、腫瘤分期等臨床信息，采用多因素Cox 回歸模型分析風(fēng)險(xiǎn)評(píng)分是否為判斷CRC不良預(yù)后的獨(dú)立因素。

1.3.4 在內(nèi)、外部驗(yàn)證集中評(píng)估預(yù)后風(fēng)險(xiǎn)評(píng)分模型的性能 采用上述計(jì)算公式對(duì)內(nèi)、外部驗(yàn)證集中每例患者的風(fēng)險(xiǎn)評(píng)分進(jìn)行計(jì)算，并據(jù)此將每個(gè)驗(yàn)證集患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。根據(jù)內(nèi)、外部驗(yàn)證集的ROC 曲線和KM 生存曲線分析結(jié)果，評(píng)估該模型在內(nèi)、外部驗(yàn)證集中的預(yù)測(cè)性能。

1.4 基因集富集分析

為驗(yàn)證GSE39582 數(shù)據(jù)集和TCGA 數(shù)據(jù)集中低風(fēng)險(xiǎn)和高風(fēng)險(xiǎn)組間的功能差異，我們利用基因集富集分析（gene set enrichment analysis，GSEA）［27］方法比較腫瘤特征基因集（hallmark gene sets）在2組間的富集度。其中，該腫瘤特征基因集來源于Molecular Signatures Database（MSigDB） 數(shù)據(jù)庫（https：//www.gsea-msigdb.org/gsea/msigdb/index.jsp），包含炎癥反應(yīng)、低氧等50 組基因集。分析輸出結(jié)果為標(biāo)準(zhǔn)化富集分?jǐn)?shù)（normalized enrichment score，NES），依據(jù)NES 的高低衡量基因集在高、低風(fēng)險(xiǎn)組間的富集程度，從而揭示高、低風(fēng)險(xiǎn)組間差異顯著的腫瘤特征信號(hào)通路或生物過程。即在P＜0.05 前提下，NES＞1 代表基因集在高風(fēng)險(xiǎn)組中顯著富集，且NES 越大富集程度越高；NES＜-1 代表基因集在低風(fēng)險(xiǎn)組中顯著富集，且|NES|越大富集程度越高。

1.5 預(yù)后風(fēng)險(xiǎn)評(píng)分模型評(píng)估免疫治療效果的應(yīng)用分析

本研究依據(jù)相關(guān)參考文獻(xiàn)，下載2 個(gè)包含生存時(shí)間、生存狀態(tài)和免疫治療效果等臨床信息的其他癌癥的數(shù)據(jù)集，分別為348 例接受程序性死亡配體-1（programmed death ligand-1，PD-L1）免疫治療的轉(zhuǎn)移性尿路上皮癌患者的基因表達(dá)數(shù)據(jù)和臨床信息［28］及49 例接受PD-1 免疫治療的黑色素瘤患者的基因表達(dá)數(shù)據(jù)和臨床信息［29］，探究CRC預(yù)后風(fēng)險(xiǎn)評(píng)分模型在免疫治療效果評(píng)估中的應(yīng)用：①分析該模型在免疫數(shù)據(jù)集中是否成立。根據(jù)“1.3.2”中的計(jì)算公式對(duì)該2個(gè)數(shù)據(jù)集中的每例患者的風(fēng)險(xiǎn)評(píng)分進(jìn)行計(jì)算，并據(jù)此將該2個(gè)數(shù)據(jù)集分別分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組。通過KM 生存曲線對(duì)2 組患者的生存率進(jìn)行分析。②采用χ2檢驗(yàn)統(tǒng)計(jì)2個(gè)數(shù)據(jù)集中不同組患者免疫治療效果間的差異，及其與風(fēng)險(xiǎn)評(píng)分之間的關(guān)系。

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

采用R 軟件（3.6.2 版本）對(duì)所有數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析。定性資料以頻數(shù)（百分比）表示，采用χ2檢驗(yàn)進(jìn)行比較。采用對(duì)數(shù)秩檢驗(yàn)（Log-Rank 法）比較所有數(shù)據(jù)集中高、低風(fēng)險(xiǎn)組患者KM 生存分析間的差異。P＜0.05 表示差異具有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 CRC數(shù)據(jù)統(tǒng)計(jì)及分析

經(jīng)整合，RNA-seq 數(shù)據(jù)集的正常結(jié)直腸組織樣本有349 例、CRC 組織樣本有471 例。來自GEO 數(shù)據(jù)庫的10個(gè)CRC數(shù)據(jù)集的基本信息見表1，因前8個(gè)數(shù)據(jù)集同時(shí)包含了CRC 組織和正常結(jié)直腸組織的樣本，故用于后續(xù)DEGs分析。

表1 GEO數(shù)據(jù)庫的10個(gè)CRC數(shù)據(jù)集的基本信息Tab 1 Basic information of ten CRC datasets from GEO database

2.2 DEGs篩選及GO功能分析

經(jīng)篩選，GEO 數(shù)據(jù)庫中前8 個(gè)CRC 數(shù)據(jù)集的DEGs數(shù)量如圖1A 所示；經(jīng)第2 次篩選后8 個(gè)CRC 數(shù)據(jù)集共有的DEGs 為962 個(gè)（上調(diào)427 個(gè)，下調(diào)535 個(gè)），其中最為顯著的10 個(gè)上調(diào)、10 個(gè)下調(diào)的DEGs 如圖1B 所示。GO功能分析結(jié)果（圖1C）顯示，上調(diào)基因主要富集在免疫細(xì)胞遷移、趨化因子介導(dǎo)的信號(hào)通路等，下調(diào)基因主要富集在負(fù)向調(diào)節(jié)細(xì)胞遷移等。

圖1 CRC數(shù)據(jù)集中的DEGs篩選及其GO功能分析Fig 1 Screening and GO functional analysis of DEGs from the CRC datasets

RNA-seq 數(shù)據(jù)集中共篩選出1 749 個(gè)DEGs（上調(diào)800個(gè)、下調(diào)949 個(gè)）。GO 功能分析結(jié)果（圖1D）顯示，多數(shù)DEGs 富集在趨化因子介導(dǎo)的信號(hào)通路、免疫細(xì)胞遷移等。

2.3 CRC預(yù)后風(fēng)險(xiǎn)評(píng)分模型構(gòu)建及其性能評(píng)價(jià)

圖2 CRC預(yù)后風(fēng)險(xiǎn)評(píng)分模型的構(gòu)建Fig 2 Construction of CRC prognostic risk score model

采用單因素Cox 回歸模型對(duì)TCGA 和GSE39582數(shù)據(jù)集共有的1 927 個(gè)DEGs 進(jìn)行分析，得到16 個(gè)與不良預(yù)后相關(guān)的基因（圖2A）。以上述基因作為參數(shù)行LASSO回歸分析，當(dāng)logλ 達(dá)到最小值（-4.11，圖2B，即黑色虛線標(biāo)注所示）時(shí)，對(duì)應(yīng)的參數(shù)為最佳建模參數(shù)，而此時(shí)有8個(gè)最佳建模參數(shù)的LASSO 回歸系數(shù)均不為0（圖2C，黑色虛線標(biāo)注所示）；對(duì)該8 個(gè)最佳建?；颍跿IMP1（TIMP metallopeptidase inhibitor 1） 、 RGL2 （ral guanine nucleotide dissociation stimulator like 2）、 SLC39A13（solute carrier family 39 member 13）、IER5L（immediate early response 5 like）、HEYL （hes related family bHLH transcription factor with YRPW motif like）、 INHBB（inhibin subunit beta B）、DMKN（dermokine）和ELFN1（extracellular leucine rich repeat and fibronectin type Ⅲdomain containing 1）］進(jìn)行多因素Cox 回歸分析，最終得到由8 個(gè)特征基因構(gòu)成的CRC 預(yù)后風(fēng)險(xiǎn)評(píng)分模型（圖2D）。風(fēng)險(xiǎn)評(píng)分計(jì)算公式為：風(fēng)險(xiǎn)評(píng)分=TIMP1 表達(dá)值×0.495+RGL2表達(dá)值×0.946+SLC39A13表達(dá)值×（-1.517）+IER5L 表達(dá)值×0.418+HEYL 表達(dá)值×0.729+INHBB 表達(dá)值×0.315+DMKN表達(dá)值×0.218+ELFN1表達(dá)值×1.226。

通過上述計(jì)算公式，獲得訓(xùn)練集中每例CRC 患者的風(fēng)險(xiǎn)評(píng)分，并據(jù)此將患者分為高風(fēng)險(xiǎn)組（N=32）和低風(fēng)險(xiǎn)組（N=187）；不同組別的患者評(píng)分分布如圖3A 所示，紅色虛線表示分組的分界線。與高風(fēng)險(xiǎn)組相比，低風(fēng)險(xiǎn)組患者的不良生存狀態(tài)較少（圖3B），且8 個(gè)建模基因的表達(dá)量更低（圖3C）。訓(xùn)練集的AUC 均值大于0.750（圖3D）；KM 分析結(jié)果顯示，與低風(fēng)險(xiǎn)組相比，高風(fēng)險(xiǎn)組患者的總體生存率更低（圖3E）。多因素Cox回歸分析的結(jié)果顯示，風(fēng)險(xiǎn)評(píng)分、年齡、腫瘤分期可以作為CRC 患者預(yù)后的獨(dú)立因素（圖3F）。以上結(jié)果證明，該模型在訓(xùn)練集中對(duì)CRC預(yù)后具有較高的預(yù)測(cè)價(jià)值。

圖3 CRC預(yù)后風(fēng)險(xiǎn)評(píng)分模型在訓(xùn)練集中的性能評(píng)估Fig 3 Performance evaluation of CRC prognostic risk score model in the training set

隨后，在內(nèi)、外部驗(yàn)證集中評(píng)估該模型的性能。內(nèi)、外部驗(yàn)證集的AUC 均值大于0.600（圖4A～D）；且該4 個(gè)驗(yàn)證集的KM生存分析顯示，與低風(fēng)險(xiǎn)組相比，高風(fēng)險(xiǎn)組患者的總體生存率更低（圖4E～H）。以上結(jié)果表明，該模型在內(nèi)、外部驗(yàn)證集中的預(yù)測(cè)能力具有中等準(zhǔn)確度。

圖4 CRC預(yù)后風(fēng)險(xiǎn)評(píng)分模型在內(nèi)、外部獨(dú)立驗(yàn)證集中的性能評(píng)估Fig 4 Performance evaluation of CRC prognostic risk score model in internal/external validation sets

2.4 高、低風(fēng)險(xiǎn)組腫瘤特征基因集的相關(guān)通路富集分析

GSEA 結(jié)果顯示，在TCGA 數(shù)據(jù)集（圖5A） 和GSE39582 數(shù)據(jù)集（圖5B）中，上皮細(xì)胞-間充質(zhì)轉(zhuǎn)化（epithelial-mesenchymal transition，EMT）生物過程、轉(zhuǎn)化生長因子β（transforming growth factor-β，TGF-β）信號(hào)通路、Kirsten 大鼠肉瘤病毒原癌基因同源產(chǎn)物（Kirsten rat sarcoma viral oncogene homolog，KRAS） 信號(hào)通路、炎癥反應(yīng)等在高風(fēng)險(xiǎn)組中具有較高的NES，而蛋白分泌、氧化磷酸化等生物過程在低風(fēng)險(xiǎn)組中具有較高的NES；從而表明，高、低風(fēng)險(xiǎn)組間的腫瘤特征基因集相關(guān)信號(hào)通路或生物過程在2個(gè)數(shù)據(jù)集中的富集均具有顯著差異。

圖5 TCGA數(shù)據(jù)集（A）和GSE39582數(shù)據(jù)集（B）中高、低風(fēng)險(xiǎn)組腫瘤特征基因集相關(guān)信號(hào)通路或生物過程的NESFig 5 NES of cancer hallmark signaling pathway or biological process between the high risk group and low risk group in the TCGA dataset(A)and GSE39582 dataset(B)

2.5 預(yù)后風(fēng)險(xiǎn)評(píng)分模型的應(yīng)用

通過KM 生存曲線進(jìn)行分別對(duì)2個(gè)數(shù)據(jù)集中高、低風(fēng)險(xiǎn)組患者進(jìn)行分析，結(jié)果（圖6）顯示，與低風(fēng)險(xiǎn)組相比，高風(fēng)險(xiǎn)組患者的生存率均較低；繼而說明，在免疫治療數(shù)據(jù)集中風(fēng)險(xiǎn)評(píng)分的高低會(huì)影響患者最終的生存狀態(tài)。因此，該CRC 預(yù)后模型在接受免疫治療的其他癌癥患者中也具有較好的評(píng)估效果。

圖6 Anti-PD-L1數(shù)據(jù)集（A）和anti-PD-1數(shù)據(jù)集（B）中高、低風(fēng)險(xiǎn)組患者的生存分析Fig 6 KM survival analyses of the patients in the high and low risk group treated with anti-PD-L1(A)and anti-PD-1(B)

為了進(jìn)一步分析患者的風(fēng)險(xiǎn)評(píng)分與免疫治療效果間的相關(guān)性，排除臨床信息中沒有免疫治療效果的患者后，分別對(duì)2個(gè)數(shù)據(jù)集中剩余的高、低風(fēng)險(xiǎn)組患者的免疫治療效果進(jìn)行分析，結(jié)果（表2）顯示低風(fēng)險(xiǎn)組患者中免疫治療效果為部分緩解（partial response，PR）和完全緩解（complete response，CR）的患者之和的占比較高，且經(jīng)χ2檢驗(yàn)分析顯示高、低風(fēng)險(xiǎn)組患者的上述免疫治療效果間差異均具有統(tǒng)計(jì)學(xué)意義（anti-PD-L1 數(shù)據(jù)集：P=0.007；anti-PD-1數(shù)據(jù)集：P=0.047）。

表2 2個(gè)數(shù)據(jù)集中高、低風(fēng)險(xiǎn)組患者的不同免疫治療效果分析Tab 2 Analysis of immunotherapy effect between the high risk group and low risk group in the two datasets

3 討論

由于缺乏適當(dāng)?shù)脑\斷方法，CRC 患者往往在晚期才能得到確診，因此針對(duì)該疾病的臨床治療及預(yù)后研究一直是學(xué)者們關(guān)注的焦點(diǎn)［30］。目前，部分研究已就CRC 患者預(yù)后模型的構(gòu)建進(jìn)行報(bào)道，如Chen等［31］構(gòu)建包含9個(gè)特征基因的CRC 預(yù)后模型，Zuo 等［32］構(gòu)建6 個(gè)特征基因的CRC 預(yù)后模型，Pagès 等［33］報(bào)道含3 個(gè)特征基因的CRC預(yù)后模型，Zhao等［34］構(gòu)建包含9個(gè)特征基因的CRC預(yù)后模型等。與之相比，本研究的CRC 預(yù)后風(fēng)險(xiǎn)評(píng)分模型存在如下優(yōu)勢(shì)：①本研究的模型構(gòu)建綜合了來自不同數(shù)據(jù)庫的轉(zhuǎn)錄組數(shù)據(jù)。相比Chen 等［31］使用的3 個(gè)GEO數(shù)據(jù)集［GSE32323（癌和癌旁樣本各17 例）、GSE74602（癌和癌旁樣本各30 例）、GSE113513（癌和癌旁樣本各14 例）］，本研究使用了更為充足的GEO 樣本量；相比Zuo 等［32］基于TCGA 數(shù)據(jù)集（647 例CRC 和51 例正常結(jié)直腸樣本） 構(gòu)建的模型，本研究不僅綜合了GEO 和TCGA 的數(shù)據(jù)，還整合了正常組織樣本數(shù)據(jù)庫GTEx 中數(shù)據(jù)使癌（471 例）和癌旁（349 例）樣本數(shù)量更為均衡。②本研究的訓(xùn)練集AUC 高于以往模型的訓(xùn)練集AUC。Chen 等［31］得出的訓(xùn)練集5 年AUC 為0.741，Pagès 等［32］分 析 的 訓(xùn) 練 集3 年AUC 為0.711、5 年AUC 為0.683，Zhao 等［34］獲 得 的 訓(xùn) 練 集3～6 年 的AUC 分 別 為0.627、0.632、0.630、0.626，上述AUC 值均低于本研究訓(xùn)練集AUC 值（1 年為0.780、2 年為0.788、3 年為0.754）。同時(shí)，本研究通過3個(gè)外部驗(yàn)證集對(duì)預(yù)后模型的預(yù)測(cè)效能進(jìn)行檢驗(yàn)（AUC 均值＞0.600），并證明該模型在其他獨(dú)立數(shù)據(jù)集中也具有中等程度的準(zhǔn)確性，而上述的前人研究并未在多個(gè)數(shù)據(jù)集中加以驗(yàn)證。③本研究對(duì)預(yù)后模型在預(yù)測(cè)癌癥患者免疫治療效果方面的應(yīng)用做了初步探索，這是其他CRC預(yù)后模型未涉足的領(lǐng)域。

構(gòu)建可靠的預(yù)后模型離不開對(duì)建模參數(shù)的逐步篩選，本研究整合了8 個(gè)GEO 數(shù)據(jù)集，并獲得了其共有的DEGs。在顯著上調(diào)的共有DEGs中，F(xiàn)OXQ1可在CRC 中促進(jìn)腫瘤相關(guān)巨噬細(xì)胞的募集引發(fā)EMT［35］；顯著下調(diào)的共有DEGs 中，ABCG2 可減輕氧化應(yīng)激和炎癥反應(yīng)，在CRC 中發(fā)揮潛在的保護(hù)作用［36］。GO 功能分析顯示，微陣列和RNA-seq數(shù)據(jù)集篩選的DEGs均富集在免疫細(xì)胞遷移、趨化因子介導(dǎo)的信號(hào)通路等。而上述富集的生物過程或信號(hào)通路與已有文獻(xiàn)研究相一致。如Waugh等［37］發(fā)現(xiàn)，在CRC 細(xì)胞中促炎趨化因子的表達(dá)較高；同時(shí)，Ackermann等［38］也證實(shí)，促炎趨化因子可促進(jìn)CRC 中嗜中性粒細(xì)胞的遷移，導(dǎo)致免疫細(xì)胞浸潤增加，且這些過程對(duì)腫瘤微環(huán)境的改變具有十分重要的作用。

本研究通過DEGs 篩選、單因素Cox 回歸分析、LASSO 回歸和多因素Cox 回歸分析，得到8 個(gè)用于建模的特征基因，即TIMP1、RGL2、SLC39A13、IER5L、HEYL、INHBB、DMKN、ELFN1。研究［39］表明，敲低TIMP1 可抑制結(jié)腸癌細(xì)胞的增殖、遷移和侵襲，并抑制CRC 中的腫瘤發(fā)生和轉(zhuǎn)移。Vigil 等［40］發(fā)現(xiàn)，RGL2 的異常過表達(dá)會(huì)促進(jìn)胰腺癌的生長。文獻(xiàn)［41］證實(shí)，與正常乳腺組織相比，乳腺癌組織中的SLC39A13 異常高表達(dá)，且與不良的生存狀態(tài)顯著相關(guān)。Ruan 等［42］發(fā)現(xiàn)，IER5L 在CRC 細(xì)胞中的表達(dá)較低。Liu 等［43］報(bào)道，敲低HEYL 可顯著減少胃癌細(xì)胞的增殖和遷移。研究［44］發(fā)現(xiàn)，CRC組織中INHBB 表達(dá)可發(fā)生異常改變。Morris 等［45］對(duì)CRC患者的結(jié)腸黏膜組織進(jìn)行測(cè)序，結(jié)果顯示相比正常的結(jié)直腸組織，CRC 患者的結(jié)腸黏膜組織中DMKN mRNA 的非翻譯區(qū)發(fā)生了異常選擇性剪接和聚腺苷酸化修飾，并推測(cè)該基因可能是CRC 的新型生物標(biāo)志物。Lei 等［46］報(bào)道ELFN1 可在CRC 細(xì)胞中發(fā)揮促增殖、抗凋亡和促遷移的功能。綜合上述研究結(jié)果，我們發(fā)現(xiàn)本研究建立的模型可能與腫瘤的發(fā)生和遷移密切相關(guān)，且這些建?；蚩赡苁悄[瘤免疫微環(huán)境中發(fā)生異常改變的關(guān)鍵環(huán)節(jié)。

同時(shí)，本研究GSEA 結(jié)果顯示高風(fēng)險(xiǎn)組患者在EMT生物過程、TGF-β 和KRAS 信號(hào)通路中具有較高的富集評(píng)分。有文獻(xiàn)報(bào)道，EMT 通路與PD-L1 之間具有雙向調(diào)節(jié)的作用［47］，且該通路和PD-1/PD-L1聯(lián)合靶向治療也是一種新的免疫治療策略［48］；另有研究發(fā)現(xiàn)，TGF-β［49］和KRAS信號(hào)通路［50］均參與了EMT的調(diào)節(jié)。綜合上述結(jié)果我們推測(cè)，高風(fēng)險(xiǎn)組患者富集評(píng)分較高的腫瘤特征通路之間可能存在相互調(diào)節(jié)，而參與這些通路的基因或許可作為高風(fēng)險(xiǎn)組患者的生物標(biāo)志物。

在預(yù)后模型的應(yīng)用領(lǐng)域，F(xiàn)riedlander 等［51］基于210例未接受過任何免疫治療的黑色素瘤樣本構(gòu)建了包括TIMP1 在內(nèi)的15 個(gè)特征基因的模型，在一定程度上可預(yù)測(cè)黑色素瘤患者的CTLA-4免疫治療效果。本研究構(gòu)建了8個(gè)特征基因的預(yù)后風(fēng)險(xiǎn)評(píng)分模型，根據(jù)風(fēng)險(xiǎn)評(píng)分將尿路上皮癌和黑色素瘤的免疫治療患者分為高、低風(fēng)險(xiǎn)組，結(jié)果顯示，高風(fēng)險(xiǎn)組患者的總體生存較差，而低風(fēng)險(xiǎn)組患者接受免疫治療的效果更好，表現(xiàn)為免疫治療效果為PR、CR的患者數(shù)量之和占比更多；上述結(jié)果表明，本研究構(gòu)建的CRC 預(yù)后風(fēng)險(xiǎn)評(píng)分模型也適用于其他癌癥，風(fēng)險(xiǎn)評(píng)分不同的患者的免疫治療效果間存在較大差異，而這種差異可能影響癌癥患者最終的生存率。因此，未來或可將本研究的預(yù)后風(fēng)險(xiǎn)評(píng)分模型用于癌癥患者的免疫效果的預(yù)測(cè)，以節(jié)約高風(fēng)險(xiǎn)患者的時(shí)間成本，優(yōu)化臨床治療方案。

綜上所述，本研究成功構(gòu)建了包含8 個(gè)特征基因的CRC 預(yù)后模型，對(duì)CRC 患者的預(yù)后具有一定的預(yù)測(cè)能力，可為挖掘CRC 免疫治療新靶點(diǎn)、優(yōu)化腫瘤免疫治療方案提供一定的參考。