高晨 吳林玉 孔寧 婁新璟 郭勇 許茂盛

肺癌是世界上發(fā)病率第二及死亡率最高的惡性腫瘤[1]。據(jù)研究報(bào)道，中國肺癌患者的5年總生存期僅為19.8%[2]。盡管手術(shù)方式、靶向治療等技術(shù)的發(fā)展明顯改善了肺癌患者的預(yù)后，但是仍然達(dá)不到預(yù)期效果[3-4]。肺腺癌是肺癌最常見的病理亞型，目前判斷肺腺癌預(yù)后的主要依據(jù)是腫瘤的TNM分期，但這種方法存在較大的差異，精確性也有待提高[5-6]。多種分子機(jī)制參與了肺腺癌的發(fā)生、發(fā)展，只有對(duì)這些機(jī)制進(jìn)行更深入的研究，才能更好地預(yù)測患者的預(yù)后，指導(dǎo)臨床治療，提高患者總生存期[7-8]。因此，臨床上迫切需要尋找更為精準(zhǔn)預(yù)測肺腺癌患者預(yù)后的分子生物標(biāo)志物。本研究擬通過基于單樣本基因富集分析（single sample gene set enrichment analysis,ssGSEA）的方法對(duì)來自癌癥基因組圖譜（The Cancer Genome Atlas,TCGA）數(shù)據(jù)庫及基因表達(dá)數(shù)據(jù)庫（Gene Expression Omnibus,GEO）中肺腺癌患者轉(zhuǎn)錄組的數(shù)據(jù)進(jìn)行系統(tǒng)性生物信息學(xué)分析并構(gòu)建肺腺癌的風(fēng)險(xiǎn)預(yù)測模型，為臨床醫(yī)師判斷肺腺癌患者總生存期提供輔助工具，同時(shí)為尋找潛在的靶向治療藥物提供參考依據(jù)。

1 資料和方法

1.1 數(shù)據(jù)獲取及整理 從基因組數(shù)據(jù)共享（Genomic Data Commons,GDC）官方網(wǎng)站（https://portal.gdc.cancer.gov）以及GEO官方網(wǎng)站（https://www.ncbi.nlm.nih.gov/geo）獲取TCGA數(shù)據(jù)庫及GSE26939數(shù)據(jù)集肺腺癌患者的轉(zhuǎn)錄組數(shù)據(jù)以及相關(guān)患者臨床數(shù)據(jù)。TCGA肺腺癌數(shù)據(jù)庫有515例患者的594個(gè)樣本，其中59個(gè)為正常樣本/癌旁樣本，535個(gè)為腫瘤樣本。GEO數(shù)據(jù)集作為外部驗(yàn)證集，數(shù)據(jù)集中有116例肺腺癌患者的116個(gè)腫瘤樣本。從免疫學(xué)數(shù)據(jù)庫和分析門戶（The Immunology Database and Analysis Portal,ImmPort）網(wǎng)站（https://www.immport.org/home）共下載免疫相關(guān)的基因1 793個(gè)。從順反組數(shù)據(jù)瀏覽器（Cistrome Data Browser,Cistrome DB）轉(zhuǎn)錄因子網(wǎng)站（http://www.cistrome.org）共下載轉(zhuǎn)錄因子318個(gè)。

1.2 ssGSEA及免疫分析 對(duì)來自TCGA數(shù)據(jù)庫肺腺癌患者的轉(zhuǎn)錄組數(shù)據(jù)，運(yùn)用ssGSEA分析及聚類分析，從而獲得樣本的免疫評(píng)分及不同免疫評(píng)分組；運(yùn)用ESTIMATE分析，獲得腫瘤樣本的腫瘤微環(huán)境分?jǐn)?shù)，并進(jìn)一步分析腫瘤微環(huán)境分?jǐn)?shù)在不同免疫評(píng)分組中的差異。運(yùn)用CIBERSORT算法計(jì)算樣本中22種腫瘤浸潤免疫細(xì)胞的比例，并分析不同免疫評(píng)分組中腫瘤浸潤免疫細(xì)胞比例的差異。在GSEA軟件（版本4.1.0）上基于不同免疫評(píng)分組進(jìn)行GSEA富集分析，獲得顯著的富集通路，并分析患者腫瘤樣本的基因表達(dá)量在不同免疫評(píng)分組中的差異，獲得在不同免疫分型組中的差異表達(dá)基因。通過取差異表達(dá)基因列表與免疫相關(guān)基因列表的交集，獲得存在差異的免疫基因。

1.3 預(yù)后基因的篩選 在得到TCGA與GEO數(shù)據(jù)集中差異免疫基因表達(dá)量后，對(duì)其進(jìn)行批次矯正。同時(shí)對(duì)患者進(jìn)行進(jìn)一步篩選，納入標(biāo)準(zhǔn)：（1）病理檢查確定為肺腺癌；（2）具有可獲取的轉(zhuǎn)錄組數(shù)據(jù)；（3）具有可獲取的臨床數(shù)據(jù)。排除標(biāo)準(zhǔn)：（1）未同時(shí)具有可獲取的轉(zhuǎn)錄組數(shù)據(jù)及臨床數(shù)據(jù)；（2）缺少相關(guān)生存數(shù)據(jù)。最終確定納入下一步研究的肺腺癌患者，其中TCGA數(shù)據(jù)庫504例，GEO數(shù)據(jù)庫115例。將來自TCGA數(shù)據(jù)集中的患者以7∶3隨機(jī)分為訓(xùn)練集356例和內(nèi)部驗(yàn)證集148例。同時(shí)將GEO數(shù)據(jù)集的患者作為外部驗(yàn)證集?；谟?xùn)練集的數(shù)據(jù)，通過單因素Cox回歸分析，篩選出肺腺癌患者預(yù)后相關(guān)的差異免疫基因。將預(yù)后相關(guān)的差異免疫基因與轉(zhuǎn)錄因子進(jìn)行相關(guān)性分析，以及在Cytoscape軟件（版本3.8.2）上進(jìn)行相關(guān)蛋白互作網(wǎng)絡(luò)的分析。

1.4 預(yù)后模型的建立 通過套索算法（least absolute shrinkage and selection operator,Lasso）回歸對(duì)篩選出的肺腺癌患者預(yù)后相關(guān)差異免疫基因進(jìn)行降維，并利用多元逐步Cox回歸模型篩選出最優(yōu)的差異免疫基因集合構(gòu)建肺腺癌的風(fēng)險(xiǎn)預(yù)測模型，并獲得每個(gè)樣本的風(fēng)險(xiǎn)分?jǐn)?shù)（Riskscore）?；讷@得Riskscore的中位數(shù)，將患者分為高風(fēng)險(xiǎn)組及低風(fēng)險(xiǎn)組。運(yùn)用內(nèi)部及外部驗(yàn)證集的數(shù)據(jù)對(duì)預(yù)測模型進(jìn)行檢驗(yàn)，并采用ROC曲線及校正曲線來顯示預(yù)測模型在訓(xùn)練組和驗(yàn)證組中預(yù)測模型的效能。采用Kaplan-Meier法對(duì)訓(xùn)練集、內(nèi)部驗(yàn)證集、外部驗(yàn)證集進(jìn)行生存分析。

1.5 聯(lián)合模型及列線圖構(gòu)建 運(yùn)用單因素及多因素Cox風(fēng)險(xiǎn)回歸，將上述獲得的Riskscore與患者的臨床特征（性別、年齡、腫瘤分期）進(jìn)行獨(dú)立預(yù)后分析，獲得肺腺癌患者獨(dú)立預(yù)后因子并構(gòu)建聯(lián)合模型。采用ROC曲線、校正曲線及列線圖分析該聯(lián)合模型的效能及臨床實(shí)用性。

1.6 統(tǒng)計(jì)學(xué)處理 采用R 4.0.5統(tǒng)計(jì)軟件。運(yùn)用的R語言包有 GSVA、limma、GSEABase、sparcl、Rtsne、ggplot2、estimate、pheatmap、reshape、reshape2、ggpubr、preprocessCore、venn、sva、survival、ggalluvial、dplyr、caret、glmnet、survminer、timeROC、rms、ggExtra、tidyverse、regplot。差異分析使用Wilcoxon秩和檢驗(yàn)，相關(guān)性分析采用Pearson相關(guān)。P＜0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 免疫評(píng)分及相關(guān)差異分析 根據(jù)腫瘤樣本的ssGSEA免疫評(píng)分及聚類分析結(jié)果，將535個(gè)腫瘤樣本分為高免疫評(píng)分組313個(gè)和低免疫評(píng)分組222個(gè)。高免疫評(píng)分組、低免疫評(píng)分組與腫瘤微環(huán)境分?jǐn)?shù)的相關(guān)性見圖1a（插頁），其中高免疫評(píng)分組中的腫瘤純度較低免疫評(píng)分組低，基質(zhì)分?jǐn)?shù)、免疫分?jǐn)?shù)及ESTIMATE分?jǐn)?shù)較低免疫評(píng)分組高，見圖1b（插頁）。運(yùn)用Wilcoxon秩和檢驗(yàn)分析基于CIBERSORT算法的12種腫瘤浸潤免疫細(xì)胞的比例在高免疫評(píng)分組和低免疫評(píng)分組中的差異有統(tǒng)計(jì)學(xué)意義，見圖1c（插頁）。在高免疫評(píng)分組中顯著富集通路有“cytokine-cytokine receptor interaction”“natural killer cell mediated cytotoxicity”“cell adhesion molecules”“T cell receptor signaling pathway”和“chemokine signaling pathway”等，見圖2a（插頁）。根據(jù)差異分析，在高免疫評(píng)分組和低免疫評(píng)分組中的差異表達(dá)基因有1 447個(gè)，見圖2b（插頁）。將差異表達(dá)基因與從ImmPort網(wǎng)站獲取的免疫相關(guān)基因取交集后，獲得免疫相關(guān)的差異表達(dá)基因有382個(gè)，見圖2c（插頁）。

圖1 肺腺癌患者免疫分析及不同免疫評(píng)分組的差異分析（a：腫瘤純度、ESTIMATE分?jǐn)?shù)、免疫分?jǐn)?shù)及基質(zhì)分?jǐn)?shù)與高免疫評(píng)分組和低免疫評(píng)分組之間的相關(guān)性，高免疫評(píng)分組中腫瘤純度較低免疫評(píng)分組低，基質(zhì)分?jǐn)?shù)、免疫分?jǐn)?shù)及ESTIMATE分?jǐn)?shù)較低免疫評(píng)分組高；b：基質(zhì)分?jǐn)?shù)、免疫分?jǐn)?shù)及ESTIMATE分?jǐn)?shù)在高免疫評(píng)分組和低免疫評(píng)分組間差異有統(tǒng)計(jì)學(xué)意義，高免疫評(píng)分組的3個(gè)分?jǐn)?shù)均比低免疫評(píng)分組高；c：基于CIBERSORT算法的12種腫瘤浸潤免疫細(xì)胞的比例在高免疫評(píng)分組和低免疫評(píng)分組中差異有統(tǒng)計(jì)學(xué)意義）

圖2 GSEA富集分析及差異免疫基因確定（a：高免疫評(píng)分組中前5個(gè)顯著富集通路有“cytokine-cytokine receptor interaction”“natural killer cell mediated cytotoxicity”“cell adhesion molecules”“ T cell receptor signaling pathway”和“chemokine signaling pathway”；b：火山圖顯示了在高免疫評(píng)分組和低免疫評(píng)分組中差異表達(dá)基因，綠色的是在高免疫評(píng)分組中低表達(dá)的基因，紅色的是在高免疫評(píng)分組中高表達(dá)的基因；c：韋恩圖顯示了綠色的是1 447個(gè)差異表達(dá)基因，粉紅色的免疫基因是1 793個(gè)，交集后免疫相關(guān)的差異表達(dá)基因有382個(gè)）

2.2 確定預(yù)后相關(guān)的差異免疫基因及其相關(guān)性分析 在對(duì)TCGA數(shù)據(jù)集及GEO數(shù)據(jù)集的轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行批次矯正且將其與上述免疫相關(guān)的差異表達(dá)基因列表取交集后，共得到219個(gè)差異表達(dá)基因的轉(zhuǎn)錄組數(shù)據(jù)。經(jīng)過進(jìn)一步篩選，共納入619例肺腺癌患者，男287例，女332例，年齡33～90（65.04±10.20）歲，患者基線資料見表1。在訓(xùn)練集中，采用單因素Cox風(fēng)險(xiǎn)回歸模型對(duì)219個(gè)基因的轉(zhuǎn)錄組數(shù)據(jù)及生存信息進(jìn)行分析后，得到31個(gè)預(yù)后相關(guān)的差異免疫基因，見圖3a（插頁），這些預(yù)后相關(guān)的差異免疫基因與轉(zhuǎn)錄因子的相關(guān)性以及蛋白互作網(wǎng)絡(luò)關(guān)系見圖3b、c（插頁）。

表1 TCGA及GEO數(shù)據(jù)集中肺腺癌患者臨床資料

圖3 預(yù)后相關(guān)的差異免疫基因及與轉(zhuǎn)錄因子相關(guān)性分析[a：森林圖顯示基于單因素Cox風(fēng)險(xiǎn)回歸分析獲得的31個(gè)預(yù)后相關(guān)的差異免疫基因；b和c：?；鶊D（b）顯示與轉(zhuǎn)錄因子具有顯著相關(guān)性的16個(gè)預(yù)后相關(guān)差異免疫基因（相關(guān)系數(shù)≥0.5,P＜0.01），PPI圖（c）顯示了其蛋白互作的關(guān)系]

2.3 預(yù)后模型的構(gòu)建及生存分析 經(jīng)過Lasso回歸及多元逐步Cox風(fēng)險(xiǎn)回歸分析降維后，得到最優(yōu)的差異免疫基因數(shù)據(jù)集并建立風(fēng)險(xiǎn)預(yù)測模型以及得到Riskscore，并根據(jù)訓(xùn)練集數(shù)據(jù)Riskscore的中位數(shù)，將訓(xùn)練集、內(nèi)部驗(yàn)證集及外部驗(yàn)證集的患者分為高風(fēng)險(xiǎn)組及低風(fēng)險(xiǎn)組，其中訓(xùn)練集高風(fēng)險(xiǎn)組178例，低風(fēng)險(xiǎn)組178例；內(nèi)部驗(yàn)證集高風(fēng)險(xiǎn)組76例，低風(fēng)險(xiǎn)組72例；外部驗(yàn)證集高風(fēng)險(xiǎn)組50例，低風(fēng)險(xiǎn)組65例。該模型的公式如下：Riskscore=EXP[（-0.211 434 889）×CX3CR1+（0.293 765 44）×IL-32+（-0.071 165 091）×SFTPD+（-0.333 423 936）×CXCR6+0.419 839 844×TAP2+（-0.269 368 749）×HLA-DOB+（-0.374 908 714）×ARG2+0.178 859 695×FURIN+（-1.040 515 441 267 49）]。該風(fēng)險(xiǎn)預(yù)測模型在對(duì)訓(xùn)練集及兩個(gè)驗(yàn)證集患者的生存預(yù)測上均具有較好的表現(xiàn)，預(yù)測患者5年總生存期的AUC分別為0.703、0.713、0.750，見圖4a-c（插頁）。同時(shí)校正曲線分析顯示該模型預(yù)測的1、3、5年患者總生存期都與實(shí)際總生存期較為一致，見圖4d-f（插頁）。此外，高風(fēng)險(xiǎn)組和低分險(xiǎn)組患者的總生存期在3個(gè)數(shù)據(jù)集中差異均有統(tǒng)計(jì)學(xué)意義，見圖4g-i（插頁）。

圖4 預(yù)測模型的ROC曲線、校正曲線分析及生存分析[a-c：訓(xùn)練集、內(nèi)部驗(yàn)證集、外部驗(yàn)證集的預(yù)測模型ROC曲線；d-f：訓(xùn)練集、內(nèi)部驗(yàn)證集、外部驗(yàn)證集的校正曲線分析，結(jié)果顯示預(yù)測結(jié)果與實(shí)際結(jié)果較為一致；g-i：在訓(xùn)練集、內(nèi)部驗(yàn)證集、外部驗(yàn)證集中，高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組患者生存曲線差異均有統(tǒng)計(jì)學(xué)意義（均P＜0.01），低風(fēng)險(xiǎn)組的患者總生存期較高風(fēng)險(xiǎn)組高]

2.4 獨(dú)立預(yù)后分析及列線圖的構(gòu)建 經(jīng)單因素及多因素Cox風(fēng)險(xiǎn)回歸的分析結(jié)果顯示，肺腺癌的分期及預(yù)測模型的Riskscore可作為肺腺癌患者的兩個(gè)獨(dú)立預(yù)后因子，見表2。肺腺癌的分期聯(lián)合預(yù)測模型的Riskscore構(gòu)建聯(lián)合模型的1、3、5年總生存期的AUC分別為0.789、0.763、0.746。校正曲線分析顯示該模型預(yù)測的1、3、5年患者總生存期與實(shí)際總生存期較為一致。同時(shí)該聯(lián)合模型臨床應(yīng)用的列線圖見圖5。

表2 單因素及多因素Cox回歸分析

圖5 聯(lián)合模型的列線圖的構(gòu)建（點(diǎn)表示1例患者，其Riskscore與分期Ⅳ期的總得分為76.9分，總生存期＜1年的概率是0.196，＜3年的概率是0.578，＜5年的概率0.84）

3 討論

ssGSEA是能對(duì)單個(gè)樣本的通路富集情況進(jìn)行量化并評(píng)分的分析方法，例如免疫相關(guān)通路等[9]。同時(shí)已有不少研究報(bào)道了腫瘤微環(huán)境中免疫細(xì)胞浸潤水平與腫瘤的發(fā)生、發(fā)展以及患者預(yù)后均相關(guān)[10-11]。而基于ssGSEA免疫分型來分析免疫相關(guān)差異表達(dá)基因在肺腺癌預(yù)后方面價(jià)值的研究尚缺乏。

本研究基于TCGA數(shù)據(jù)集中肺腺癌患者樣本基因表達(dá)數(shù)據(jù)的ssGSEA免疫評(píng)分及Cox風(fēng)險(xiǎn)回歸分析，獲得肺腺癌預(yù)后相關(guān)的免疫差異表達(dá)基因。然后進(jìn)行Lasso回歸及多元逐步Cox風(fēng)險(xiǎn)回歸分析降維，構(gòu)建肺腺癌預(yù)后預(yù)測模型獲得Riskscore，并用GEO數(shù)據(jù)集進(jìn)行驗(yàn)證。利用Riskscore預(yù)測患者5年總生存期的AUC值均＞0.7。將Riskscore與臨床特征進(jìn)行獨(dú)立預(yù)后因子分析，并建立列線圖，列線圖1、3、5年的AUC也均＞0.7。上述結(jié)果表明Riskscore及列線圖在預(yù)測肺腺癌總生存期上都具有良好的效能。

本研究結(jié)果顯示肺腺癌高免疫評(píng)分組的富集通路主要在“cytokine-cytokine receptor interaction”“natural killer cell mediated cytotoxicity”和“cell adhesion molecules”。而國內(nèi)外也已有文獻(xiàn)報(bào)道這些通路在非小細(xì)胞肺癌中的發(fā)生、發(fā)展中起重要作用[12-14]。例如，Zheng等[12]研究結(jié)果表明在高免疫浸潤組的肺腺癌腫瘤樣本也主要富集在“cytokine-cytokine receptor interaction”。Zhang等[13]研究也發(fā)現(xiàn)EPHA5突變會(huì)破壞自然殺傷細(xì)胞介導(dǎo)的細(xì)胞毒性對(duì)非小細(xì)胞肺癌的作用，促進(jìn)癌細(xì)胞遷移。Li等[14]研究表明高免疫評(píng)分組和低免疫評(píng)分組的差異基因富集通路也包括“cytokine-cytokine receptor interaction”和“cell adhesion molecules”。上述研究結(jié)果證實(shí)了這些通路可以作為肺腺癌潛在的免疫治療新靶點(diǎn)，進(jìn)而提高肺腺癌患者預(yù)后水平。

本研究中，在對(duì)預(yù)后相關(guān)的免疫差異表達(dá)基因進(jìn)行Lasso回歸及多元逐步Cox風(fēng)險(xiǎn)回歸分析降維后，獲得了最優(yōu)的8個(gè)基因集合，其中權(quán)重最高的3個(gè)基因分別是TAP2、ARG2和CXCR6。這些基因在非小細(xì)胞肺癌以及肺腺癌進(jìn)展方面的作用也已經(jīng)有一些研究報(bào)道[15-17]。Liu等[15]研究結(jié)果表明TAP2的基因多態(tài)性與非小細(xì)胞肺癌的進(jìn)展存在潛在的關(guān)系。Giatromanolaki等[16]研究結(jié)果發(fā)現(xiàn)ARG2主要在癌相關(guān)成纖維細(xì)胞中表達(dá)。說明這些基因具有預(yù)測肺腺癌患者預(yù)后的能力。但是，也有研究結(jié)果發(fā)現(xiàn)癌細(xì)胞中表達(dá)的CXCR16并沒有表現(xiàn)出對(duì)非小細(xì)胞肺癌預(yù)后的影響[17]。這可能是該研究納入的是非小細(xì)胞肺癌患者，與本研究中納入的肺腺癌患者不同所致，這需要進(jìn)一步深入的研究來驗(yàn)證。

通過類似的方法，Wang等[18]研究也是基于ssGSEA算法分析出的免疫相關(guān)差異表達(dá)的長鏈非編碼RNA（lncRNA）來構(gòu)建肺腺癌預(yù)后的預(yù)測模型，其訓(xùn)練集預(yù)測5年總生存期的AUC達(dá)到0.63，驗(yàn)證集預(yù)測5年總生存期的AUC達(dá)到0.667。本研究則基于ssGSEA分析構(gòu)建肺腺癌患者風(fēng)險(xiǎn)預(yù)測模型并獲得Riskscore，其在訓(xùn)練集、內(nèi)部驗(yàn)證集及外部驗(yàn)證集5年總生存期的AUC分別為0.703、0.713、0.750。同時(shí)將該Riskscore與臨床特征進(jìn)行獨(dú)立預(yù)后分析并構(gòu)建聯(lián)合模型列線圖，列線圖5年總生存期的AUC達(dá)到0.746。Riskscore與列線圖的AUC均較高且均＞0.7，Riskscore的HR值也較臨床分期高，均進(jìn)一步表明了Riskscore在預(yù)測肺腺癌預(yù)后方面的能力較好，同時(shí)該模型也具有較好的泛化性。通過校正曲線分析也進(jìn)一步顯示了預(yù)測模型的結(jié)果與實(shí)際總生存期較為一致。

本研究存在一定的局限性。首先，研究納入的患者例數(shù)有限，僅用了來自于TCGA數(shù)據(jù)庫及GEO數(shù)據(jù)集，未來還需要多中心的數(shù)據(jù)來進(jìn)行驗(yàn)證。其次，本研究并未對(duì)篩選的風(fēng)險(xiǎn)基因進(jìn)行相關(guān)機(jī)制研究，未來需開展基礎(chǔ)實(shí)驗(yàn)進(jìn)行進(jìn)一步分析。

綜上所述，本研究基于ssGSEA分析獲得的免疫相關(guān)差異表達(dá)基因構(gòu)建了肺腺癌風(fēng)險(xiǎn)預(yù)測模型，并聯(lián)合臨床特征繪制了列線圖，以期能夠輔助臨床醫(yī)師對(duì)肺腺癌患者總生存期進(jìn)行判斷。