鄧彪 陳春源 黃裕鋒 陳心銘 羅連響 梁柱（廣東醫(yī)科大學(xué)附屬醫(yī)院心胸外科，湛江 524001）

肺癌是全世界癌癥相關(guān)死亡的主要原因之一，盡管近年肺癌病死率有所下降，但2017年導(dǎo)致的死亡病例數(shù)仍多于乳腺癌、前列腺癌、結(jié)直腸癌和顱內(nèi)腫瘤死亡病例總和［1］。肺腺癌（lung adenocarcinoma，LUAD）是非小細(xì)胞肺癌最普遍的組織類型，與其他肺癌亞型相比，LUAD與基因組變化關(guān)系密切，具有高度異質(zhì)性［2-3］。半數(shù)以上LUAD診斷時(shí)已是局部晚期或已轉(zhuǎn)移，預(yù)后較差，五年生存率較低，開(kāi)發(fā)一個(gè)更有效的LUAD預(yù)后預(yù)測(cè)工具迫在眉睫［4-5］。癌細(xì)胞新陳代謝是癌癥的標(biāo)志之一，可維持癌細(xì)胞生長(zhǎng)［6］。近年研究表明癌細(xì)胞新陳代謝改變對(duì)癌癥生長(zhǎng)和發(fā)展起重要作用，如LUAD預(yù)后與癌細(xì)胞新陳代謝密切相關(guān)［7-10］。因此，系統(tǒng)研究代謝相關(guān)基 因（metabolism related genes，MRGs）對(duì) 改 善LUAD治療和預(yù)后至關(guān)重要［11-12］。腫瘤免疫微環(huán)境（tumor immune microenvironment，TIME）中的免疫細(xì)胞通過(guò)分泌代謝決定因子維持腫瘤生長(zhǎng)，其類型和豐度對(duì)腫瘤進(jìn)展和免疫治療具有顯著影響［13-14］。但系統(tǒng)分析MRGs、LUAD預(yù)后與TIME關(guān)系的研究較少。本研究旨在建立基于MRGs的LUAD預(yù)測(cè)模型，并分析LUAD患者TIME與MRGs的關(guān)系，通過(guò)生物信息學(xué)綜合分析方法研究TCGA（The Cancer Genome Atlas）和GEO（Gene Expression Omnibus）數(shù)據(jù)庫(kù)的LUAD基因表達(dá)數(shù)據(jù)和臨床信息，建立一個(gè)性能穩(wěn)定的風(fēng)險(xiǎn)預(yù)測(cè)模型和列線圖，通過(guò)GO和KEGG功能富集分析探索與模型MRGs密切相關(guān)的通路，采用免疫評(píng)分相關(guān)性分析探索高低風(fēng)險(xiǎn)患者免疫細(xì)胞和免疫相關(guān)通路差異，探討MRGs與TIME的關(guān)系。

1 資料與方法

1.1樣本基因組表達(dá)數(shù)據(jù)和臨床信息研究流程如圖1所示。從TCGA數(shù)據(jù)庫(kù)（https：//ancergenome.nih.gov/）在線下載LUAD患者訓(xùn)練組轉(zhuǎn)錄數(shù)據(jù)和臨床數(shù)據(jù)（轉(zhuǎn)錄組數(shù)據(jù)共551個(gè)樣本，56 753個(gè)基因，其中腫瘤樣本497例，癌旁樣本54例；臨床數(shù)據(jù)共522個(gè)樣本，包含生存時(shí)間、生存狀態(tài)、年齡、性別、分期和T、N、M分期等臨床信息），從GEO數(shù)據(jù)庫(kù)（https：//www.ncbi.nlm.nih.gov/geo/）在線下載包括患者臨床信息和基因表達(dá)的芯片數(shù)據(jù)GSE31210，包括246個(gè)早期LUAD樣本的54 675個(gè)基因，其中226個(gè)樣本包含生存時(shí)間、生存狀態(tài)、年齡、性別、分期和吸煙史等在內(nèi)的完整臨床信息。納入本研究的LUAD患者臨床數(shù)據(jù)如表1所示。從分子簽名數(shù)據(jù)庫(kù)（http：//www.gsea-msigdb.org/gsea/index.jsp）下載959個(gè)MRGs［15-16］，篩選生存時(shí)間＞30 d的樣本進(jìn)行生存分析。

表1 納入本研究的LUAD患者臨床特征［例（%）］Tab.1 Clinical characteristics of LUAD patients included in this study［n（%）］

圖1 研究流程圖Fig.1 Flow chart

1.2MRGs表達(dá)矩陣提取將從TCGA下載的LUAD表達(dá)數(shù)據(jù)去除在所有樣本中均呈低表達(dá)的基因（在所有樣本的平均表達(dá)＜0）并進(jìn)行標(biāo)準(zhǔn)化處理，與從分子簽名數(shù)據(jù)庫(kù)下載的959個(gè)MRGs取交集，提取出949個(gè)交集MRGs的表達(dá)數(shù)據(jù)。

1.3篩選預(yù)后相關(guān)差異MRGs將TCGA訓(xùn)練隊(duì)列的腫瘤組織和正常組織分類，通過(guò)R語(yǔ)言“l(fā)imma”包在提取出交集MRGs的表達(dá)數(shù)據(jù)中進(jìn)行基因差異表達(dá)分析，以fdrFilter=0.05和logFCfilter=2為篩選界值，得到100個(gè)顯著差異表達(dá)MRGs，繪制差異MRGs火山圖。將TCGA訓(xùn)練隊(duì)列中結(jié)合臨床信息的LUAD患者表達(dá)數(shù)據(jù)進(jìn)行單因素COX分析得到209個(gè)預(yù)后相關(guān)MRGs，將顯著差異表達(dá)的MRGs與預(yù)后相關(guān)MRGs取交集，得到32個(gè)預(yù)后相關(guān)差異MRGs，以韋恩圖顯示交集過(guò)程并繪制交集基因在正常組織和LUAD腫瘤組織中的差異表達(dá)熱圖。將GEO隊(duì)列的基因表達(dá)信息與其臨床信息合并，得到驗(yàn)證組表達(dá)與生存數(shù)據(jù)，用于模型驗(yàn)證。所有分析采用R軟件（4.0.3）完成。采用String PPI網(wǎng)絡(luò)在線工 具（https：//string-db.org/）分 析 預(yù) 后 相 關(guān) 差 異MRGs的關(guān)系［17］，并輸出蛋白互作網(wǎng)絡(luò)圖。

1.4獲取共表達(dá)基因?qū)SE31210樣本得到的246個(gè)早期LUAD患者表達(dá)數(shù)據(jù)進(jìn)行處理，包括去除表達(dá)量無(wú)效值、轉(zhuǎn)換ID、去除低表達(dá)和去重，得到含246個(gè)LUAD樣本的15 255個(gè)基因GEO驗(yàn)證隊(duì)列表達(dá)數(shù)據(jù)。將GEO驗(yàn)證隊(duì)列基因與TCGA訓(xùn)練隊(duì)列預(yù)后相關(guān)差異MRGs取交集，得到19個(gè)兩隊(duì)列中共表達(dá)的MRGs，并分別在TCGA隊(duì)列與GEO隊(duì)列中提取共表達(dá)MRGs的表達(dá)矩陣。

1.5MRGs風(fēng)險(xiǎn)模型在TCGA共表達(dá)MRGs的表達(dá)矩陣中，經(jīng)Lasso回歸分析計(jì)算出各MRGs回歸系數(shù)β，篩選出12個(gè)建立風(fēng)險(xiǎn)預(yù)測(cè)模型的最佳候選MRGs，用候選基因表達(dá)水平與其回歸系數(shù)β乘積之和建立風(fēng)險(xiǎn)預(yù)測(cè)模型并繪制Lasso模型圖。風(fēng)險(xiǎn)值=基因1表達(dá)×β1+基因2表達(dá)×β2+…+基因n表達(dá)×βn。

1.6預(yù)后潛力分析以風(fēng)險(xiǎn)值中位數(shù)為界值，分別將TCGA和GEO隊(duì)列中LUAD患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，Kaplan-Meier（Km）曲線預(yù)測(cè)高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)患者總生存率差異，ROC曲線檢驗(yàn)風(fēng)險(xiǎn)模型的預(yù)后潛力和意義，AUC值評(píng)估模型預(yù)測(cè)效能。本研究采用R語(yǔ)言“Survive”包比較最優(yōu)基因在LUAD和正常組織中的表達(dá)差異。

1.7GO和KEGG富集分析采用R語(yǔ)言的R包“clusterProfiler”“org.Hs.eg.db”和“enrichplot”進(jìn)行GO和KEGG富集分析，探討不同分子機(jī)制及高、低風(fēng)險(xiǎn)患者差異。GO分析包括生物過(guò)程、細(xì)胞成分和分子功能。

1.8關(guān)鍵MRGs表達(dá)分析為分析構(gòu)建模型MRGs表達(dá)以進(jìn)一步驗(yàn)證上述分析結(jié)果，本研究采用HPA（The Human Protein Atlas）在線數(shù)據(jù)庫(kù)（https：//www.proteinatlas.org/）對(duì)比TCGA數(shù)據(jù)庫(kù)中癌組織和鄰近正常肺組織的病理切片［18］，驗(yàn)證這些MRGs在翻譯水平的差異表達(dá)情況。

1.9TIME分析不同LUAD患者生存差異可能由TIME的復(fù)雜性導(dǎo)致，因此本研究采用單樣本基因集富集分析（ssGSEA）算法分析一些免疫細(xì)胞和免疫相關(guān)通路在高、低風(fēng)險(xiǎn)患者中的比例。通過(guò)采用R語(yǔ)言“GSVA”包和“GSEABase”包對(duì)TCGA和GEO隊(duì)列中高、低風(fēng)險(xiǎn)樣本進(jìn)行免疫評(píng)分，用“l(fā)imma”“ggpubr”和“reshape2”包進(jìn)行免疫評(píng)分相關(guān)性分析，采用ssGSEA算法進(jìn)一步量化免疫細(xì)胞浸潤(rùn)水平和免疫相關(guān)通路活性，包括16種免疫細(xì)胞：aDCs、B細(xì)胞、CD8+T細(xì)胞、DCs、巨噬細(xì)胞、肥大細(xì)胞、中性粒細(xì)胞、自然殺傷細(xì)胞、pDCs、T輔助細(xì)胞、Tfh細(xì)胞、Th1細(xì)胞、Th2細(xì)胞、腫瘤浸潤(rùn)的淋巴細(xì)胞、調(diào)節(jié)性T細(xì)胞、iDCs，13種免疫相關(guān)活性通路：APC共抑制、APC共刺激、CCR、檢查點(diǎn)、溶細(xì)胞活性、HLA、驗(yàn)證促進(jìn)、MHCⅠ、副炎癥、T細(xì)胞共抑制、T細(xì)胞共刺激、Ⅰ型IFN反應(yīng)、Ⅱ型IFN反應(yīng)。

1.10統(tǒng)計(jì)學(xué)分析采用t檢驗(yàn)確定腫瘤組織和正常組織差異表達(dá)MRGs。采用卡方檢驗(yàn)比較各組比例構(gòu)成差異。Mann-Whitney檢驗(yàn)評(píng)估風(fēng)險(xiǎn)組間免疫細(xì)胞或途徑ssGSEA得分差異，BH法調(diào)整P值。單因素和多因素COX回歸分析確定OS的獨(dú)立預(yù)后因素。Lasso回歸分析計(jì)算候選基因回歸系數(shù)β，Km分析和對(duì)數(shù)秩和檢驗(yàn)（Log-rank test）比較各組OS，所有統(tǒng)計(jì)分析均采用R軟件（4.0.3版本）完成。GSEA分析各通路矯正富集評(píng)分（NES）與富集程度（NOM p-val），所有分析均以P＜0.05為差異具有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1MRGs在LUAD和正常組織樣本存在顯著差異從TCGA中得到所有LUAD患者mRNA表達(dá)數(shù)據(jù)和臨床資料，在分子簽名數(shù)據(jù)庫(kù)（http：//www.gsea-msigdb.org/gsea/index.jsp）提取出959個(gè)MRGs。將從TCGA下載的LUAD表達(dá)數(shù)據(jù)去除在所有樣本中均呈低表達(dá)的基因并進(jìn)行標(biāo)準(zhǔn)化處理，與959個(gè)MRGs取交集，提取出949個(gè)交集MRGs表達(dá)數(shù)據(jù)。TCGA訓(xùn)練隊(duì)列中腫瘤組織和正常組織分類后，通過(guò)R語(yǔ)言“l(fā)imma”包在提取出交集MRGs表達(dá)數(shù)據(jù)中進(jìn)行差異表達(dá)分析，以fdrFilter=0.05和logFCfilter=2作為篩選界值，得到100個(gè)在腫瘤和正常組織顯著差異表達(dá)MRGs，繪制差異MRGs火山圖（圖2A）。在TCGA訓(xùn)練隊(duì)列中結(jié)合臨床信息的LUAD患者基因進(jìn)行單因素COX分析顯示，其中209個(gè)與LUAD患者OS相關(guān)。將預(yù)后相關(guān)MRGs與顯著差異表達(dá)的MRGs取交集，得到32個(gè)與LUAD患者OS相關(guān)的差異MRGs（圖2B），韋恩圖顯示32個(gè)與預(yù)后相關(guān)的差異MRGs交集過(guò)程，并繪制了交集基因在正常組織和LUAD腫瘤組織中的差異表達(dá)熱圖（圖2C），樹(shù)狀圖（圖2D）顯示單因素COX分析32個(gè)基因的P值和風(fēng)險(xiǎn)比（HR）：16個(gè)MRGs為保護(hù)基因（HR＜1），16個(gè)MRGs為危險(xiǎn)基因（HR＞1）。

圖2 MRGs在LUAD和正常組織樣本存在顯著差異Fig.2 MRGs were significantly different between LUAD and normal tissues

2.2基于MRGs的預(yù)后模型 對(duì)TCGA-LUAD訓(xùn)練隊(duì)列中12個(gè)模型基因表達(dá)的關(guān)系進(jìn)行評(píng)估，采用String PPI網(wǎng)絡(luò)在線工具（https：//string-db.org/）繪制蛋白互作網(wǎng)絡(luò)圖（圖3A），并制作基因相關(guān)性網(wǎng)絡(luò)（圖3B），結(jié)果顯示32個(gè)MRGs存在相互作用。相關(guān)性分析表明，F(xiàn)MO2與RRM2、TK1、PFKP呈顯著負(fù)相關(guān)。將GEO隊(duì)列基因表達(dá)信息與其臨床信息合并，得到驗(yàn)證組表達(dá)與生存數(shù)據(jù)，去除表達(dá)量無(wú)效值、轉(zhuǎn)換ID、去除低表達(dá)和去重等處理后得到15 255個(gè)基因，與TCGA中32個(gè)預(yù)后相關(guān)差異MRGs取交集后得到19個(gè)共表達(dá)基因。在TCGA表達(dá)數(shù)據(jù)中提取訓(xùn)練組共表達(dá)基因表達(dá)矩陣。為整合這些關(guān)鍵基因作用，并確定對(duì)生存結(jié)果最重要的基因，在訓(xùn)練訓(xùn)練組共表達(dá)矩陣中將上述19個(gè)基因進(jìn)行Lasso回歸分析，根據(jù)Lasso回歸結(jié)果確定12個(gè)基因（TK1、GCLC、PFKP、CA9、CEL、CA4、RRM2、CYP2D6、HAL、GSTA3、FMO2和ENO3）構(gòu)建風(fēng)險(xiǎn)預(yù)測(cè)模型（圖3C、D），12個(gè)基因全名和風(fēng)險(xiǎn)系數(shù)見(jiàn)表2。模型風(fēng)險(xiǎn)值=（0.016 789 869 774 756 4×TK1表達(dá)量）+（0.011 750 257 755 730 3×GCLC表 達(dá) 量）+（0.057 003 517 487 642 8×PFKP表 達(dá) 量）+（0.033 715 712 604 503 8×CA9表 達(dá) 量）+（-0.122 404 966 481 127×CEL表 達(dá) 量）+（-0.060 890 413 642 635 5×CA4表 達(dá) 量）+（0.144 614 619 025 540 0×RRM2表 達(dá) 量）+（-0.038 521 995 672 867 4×CYP2D6表 達(dá) 量）+（0.081 884 526 436 414 1×HAL表 達(dá) 量）+（-0.311 751 816 268 740 0×GSTA3表 達(dá) 量）+（-0.054 444 453 536 976 8×FMO2表 達(dá) 量）+（-0.096 587 758 173 598 2×ENO3表達(dá)量）。

圖3 預(yù)后模型MRGs篩選Fig.3 Screening of MRGs for prognostic model

表2 預(yù)測(cè)模型中的基因Tab.2 Genes in prediction model

2.3風(fēng)險(xiǎn)評(píng)分與患者預(yù)后的關(guān)系根據(jù)轉(zhuǎn)錄矩陣中基因表達(dá)數(shù)據(jù)與臨床數(shù)據(jù)（TCGA隊(duì)列中共457個(gè)有完整臨床信息的樣本被納入風(fēng)險(xiǎn)預(yù)后分析）計(jì)算TCGA訓(xùn)練集中各LUAD患者風(fēng)險(xiǎn)分?jǐn)?shù)，將風(fēng)險(xiǎn)分?jǐn)?shù)按遞增順序排列（圖4A），以中位數(shù)將訓(xùn)練集中所有患者分為高風(fēng)險(xiǎn)組（n=228）和低風(fēng)險(xiǎn)組（n=229）。

圖4 TCGA-LUAD隊(duì)列中12基因預(yù)測(cè)模型風(fēng)險(xiǎn)評(píng)分分析Fig.4 Risk score analysis of 12-gene predictive model in TCGA-LUAD cohort

2.4風(fēng)險(xiǎn)評(píng)分是獨(dú)立預(yù)后指標(biāo) 為探討風(fēng)險(xiǎn)分?jǐn)?shù)和臨床特征與患者生存的相關(guān)性，本研究進(jìn)行了單因素和多因素COX風(fēng)險(xiǎn)回歸分析，以風(fēng)險(xiǎn)分?jǐn)?shù)、性別、年齡、臨床分期、T、N和M分期為協(xié)變量，評(píng)估以上指標(biāo)預(yù)測(cè)預(yù)后的潛力。單因素COX回歸分析（圖6A）結(jié)果顯示：臨床分期（HR：1.597；95%CI：1.361～1.873；P＜0.001）、T分期（HR：1.593；95%CI：1.304～1.947；P＜0.001）、N分期（HR：1.959；95%CI：1.101～圖4B顯示了訓(xùn)練集中457例患者風(fēng)險(xiǎn)評(píng)分、OS（以年為單位）和生存狀態(tài)，按風(fēng)險(xiǎn)評(píng)分遞增順序排列，高風(fēng)險(xiǎn)分?jǐn)?shù)患者病死率高于低風(fēng)險(xiǎn)分?jǐn)?shù)患者。對(duì)所有患者進(jìn)行主成分分析（PCA），PCA圖和t-SNE圖顯示，不同風(fēng)險(xiǎn)組患者被分為兩個(gè)方向（圖4C、D）。Km分析顯示高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組患者總體預(yù)后差異顯著（Log-ranking檢驗(yàn)，P＜0.001，圖4E），高風(fēng)險(xiǎn)組患者存活率明顯低于低風(fēng)險(xiǎn)組。12個(gè)基因差異表達(dá)熱圖（圖4F）顯示，各基因在高低風(fēng)險(xiǎn)組間表達(dá)存在差異。ROC曲線評(píng)估12個(gè)基因在LUAD患者預(yù)后預(yù)測(cè)的有效性顯示：風(fēng)險(xiǎn)評(píng)分預(yù)測(cè)三、四、五年生存率AUC分別為0.726、0.739和0.688（圖5A），對(duì)四年生存率預(yù)測(cè)方面，風(fēng)險(xiǎn)評(píng)分AUC明顯高于其他臨床特征AUC（圖5B），提示該模型預(yù)測(cè)性能良好。將本研究MRGs模型與類似模型進(jìn)行比較［19-20］，其模型三、四年ROC AUC分別為0.567和0.587（圖5D），三、四、五年ROC AUC分別為0.623、0.696和0.700（圖5E），表明本課題組開(kāi)發(fā)的模型預(yù)測(cè)性能更穩(wěn)定。3.486；P＜0.001）、M分 期（HR：1.959；95%CI：1.101～3.486；P=0.022）和風(fēng)險(xiǎn)評(píng)分（HR：3.104；95%CI：2.306～4.178；P＜0.001）與患者生存相關(guān)，多因素COX回歸分析（圖6B）結(jié)果顯示：風(fēng)險(xiǎn)評(píng)分（HR：2.676；95%CI：1.964～3.647；P＜0.001）與患者生存相關(guān)。風(fēng)險(xiǎn)評(píng)分在單因素分析和多因素分析中均具有顯著獨(dú)立預(yù)后價(jià)值（P＜0.001），證明12基因標(biāo)記對(duì)生存具有獨(dú)立預(yù)測(cè)價(jià)值。表明基于12基因標(biāo)記的風(fēng)險(xiǎn)評(píng)分對(duì)預(yù)測(cè)LUAD患者預(yù)后具有重要意義?；赥CGA-LUAD隊(duì)列進(jìn)一步構(gòu)建了臨床診療過(guò)程中實(shí)用性更強(qiáng)的列線圖（圖6C）及其校準(zhǔn)曲線（圖6D～F），納入4個(gè)預(yù)后因素（年齡、性別、分期和風(fēng)險(xiǎn)評(píng)分），根據(jù)影響生存的預(yù)后因素比例對(duì)患者進(jìn)行評(píng)分，預(yù)測(cè)患者一、三、五年生存率，列線圖預(yù)測(cè)一、三、五年生存的校準(zhǔn)曲線與對(duì)角線能較好重合，顯示出列線圖較穩(wěn)定的預(yù)測(cè)性能。

圖5 ROC曲線驗(yàn)證風(fēng)險(xiǎn)模型預(yù)測(cè)生存率的穩(wěn)定性［19-20］Fig.5 ROC curve verifies stability of risk model predicting survival rate［19-20］

圖6 獨(dú)立預(yù)后分析與列線圖Fig.6 Independent prognostic analysis and nomogram

2.5風(fēng)險(xiǎn)簽名驗(yàn)證為進(jìn)一步探討預(yù)測(cè)性能，通過(guò)GSE31210中226個(gè)樣本進(jìn)行外部驗(yàn)證，采用上述方法計(jì)算GEO驗(yàn)證集所有LUAD患者風(fēng)險(xiǎn)分?jǐn)?shù)，將風(fēng)險(xiǎn)分?jǐn)?shù)按遞增順序排列（圖7A），以中位數(shù)將訓(xùn)練集中所有患者分為高風(fēng)險(xiǎn)組（n=113）和低風(fēng)險(xiǎn)組（n=113）。圖7B顯示驗(yàn)證集中226例患者風(fēng)險(xiǎn)評(píng)分、OS（以年為單位）和生存狀態(tài)，按風(fēng)險(xiǎn)評(píng)分遞增順序排列，高風(fēng)險(xiǎn)組有更多死亡患者。所有患者PCA和t-SNE分析顯示，不同風(fēng)險(xiǎn)組患者仍可被較明顯地分為兩個(gè)方向（圖7C、D）。Km分析（圖7E）同樣顯示：高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組患者總體預(yù)后差異顯著（Log-ranking檢驗(yàn)，P＜0.05），圖7F差異熱圖顯示GEO驗(yàn)證組各基因在高低風(fēng)險(xiǎn)組間差異表達(dá)。GEO驗(yàn)證集中風(fēng)險(xiǎn)評(píng)分預(yù)測(cè)五年生存率ROC曲線AUC也達(dá)到0.681（圖5C），高于其他臨床特征AUC，驗(yàn)證了模型的可靠預(yù)后性能。

圖7 GEO-LUAD隊(duì)列中12基因預(yù)測(cè)模型風(fēng)險(xiǎn)評(píng)分分析Fig.7 Risk score analysis of 12-gene predictive model in GEO-LUAD cohort

2.6通路和功能富集分析 為闡明與風(fēng)險(xiǎn)評(píng)分相關(guān)的生物學(xué)功能和通路，基于TCGA-LUAD隊(duì)列高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組差異表達(dá)基因，對(duì)GO和KEGG途徑進(jìn)行富集分析，GO分析顯示（圖8A），生物學(xué)過(guò)程（BP）中，TCGA隊(duì)列的風(fēng)險(xiǎn)組間差異表達(dá)基因主要富集于小分子分解過(guò)程和嘌呤復(fù)合物代謝過(guò)程，細(xì)胞成分（CC）方面，差異表達(dá)基因主要富集于線粒體基質(zhì)轉(zhuǎn)移酶復(fù)合物和含磷基團(tuán)。分子功能（MF）方面，差異表達(dá)基因主要富集于磷酸酯水解酶活性和裂解酶活性。KEGG途徑分析顯示，差異基因在嘌呤代謝、輔因子生物合成及碳代謝等通路顯著富集（圖8B）。

圖8 GO和KEGG富集分析Fig.8 GO and KEGG enrichment analysis

2.7TIME為進(jìn)一步討論高風(fēng)險(xiǎn)人群和低風(fēng)險(xiǎn)人群TIME差異，評(píng)估了免疫細(xì)胞和免疫相關(guān)功能通路得分。TCGA隊(duì)列高風(fēng)險(xiǎn)組中，NK細(xì)胞、Th1細(xì)胞、Th2細(xì)胞的免疫細(xì)胞亞群顯著上調(diào)（調(diào)整后P＜0.05）；B細(xì)胞、肥大細(xì)胞、中性粒細(xì)胞、T輔助細(xì)胞、iDCs的免疫細(xì)胞亞群顯著下調(diào)（調(diào)整后P＜0.05，圖9A）。免疫相關(guān)通路方面，APC共刺激、炎癥加劇、Ⅰ型主要組織相容性復(fù)合體和副炎癥顯著上調(diào)（校正P＜0.05），人類白細(xì)胞抗原和Ⅱ型IFN反應(yīng)顯著下調(diào)（校正P＜0.05，圖9B）。通過(guò)GSE31210隊(duì)列高風(fēng)險(xiǎn)組驗(yàn)證，巨噬細(xì)胞、Th2細(xì)胞、Treg及檢查點(diǎn)、副炎癥、T細(xì)胞共刺激、Ⅰ型IFN反應(yīng)等免疫功能上調(diào)，肥大細(xì)胞、iDCs細(xì)胞和Ⅱ型IFN反應(yīng)下調(diào)（P＜0.05，圖9C、D）。

2.8MRGs在翻譯水平的差異HPA數(shù)據(jù)庫(kù)中免疫組化結(jié)果顯示：CA4在腫瘤組織中表達(dá)的蛋白被中度染色，CA9、PFKP在腫瘤組織中表達(dá)的蛋白被重度染色，RRM2在腫瘤組織中表達(dá)的蛋白也被輕度染色，而4個(gè)基因在正常肺組織中均未找到被染色的相應(yīng)蛋白（圖10）。

3 討論

近年高通量技術(shù)和生物信息學(xué)方法迅速發(fā)展，差異基因和基于特定通路的基因簽名在預(yù)測(cè)LUAD預(yù)后中的研究與日俱增，MRGs與癌癥預(yù)后關(guān)系的研究備受關(guān)注［21-25］。LUAD是非小細(xì)胞肺癌最常見(jiàn)亞型，具有顯著預(yù)后差異［11-12，26］。精準(zhǔn)預(yù)測(cè)LUAD生存率使開(kāi)發(fā)穩(wěn)定的預(yù)后指標(biāo)極為重要。本研究揭示了MRGs在LUAD中的作用及臨床意義，驗(yàn)證了12個(gè)MRGs與LUAD發(fā)生發(fā)展的密切關(guān)系，構(gòu)建了一個(gè)新的MRGs預(yù)后模型，討論了該模型作為臨床預(yù)測(cè)指標(biāo)的價(jià)值。

本研究從TCGA-LUAD數(shù)據(jù)集共檢測(cè)到100個(gè)差異表達(dá)MRGs，通過(guò)LUAD單因素COX回歸分析，發(fā)現(xiàn)了32個(gè)與OS顯著相關(guān)的MRGs。Lasso回歸分析后建立了一個(gè)由12個(gè)MRGs（TK1、GCLC、PFKP、

CA9、CEL、CA4、RRM2、CYP2D6、HAL、GSTA3、FMO2和ENO3）組成的新的預(yù)后標(biāo)志。根據(jù)預(yù)后標(biāo)志的風(fēng)險(xiǎn)評(píng)分將LUAD患者分為高風(fēng)險(xiǎn)組和低風(fēng)險(xiǎn)組，高風(fēng)險(xiǎn)組患者OS明顯低于低風(fēng)險(xiǎn)組。隨時(shí)間變化的ROC分析顯示，一、三、五年的AUC分別為0.726、0.739和0.688，提示此預(yù)后標(biāo)志具有良好的敏感性和特異性。進(jìn)一步分析表明，調(diào)整到其他臨床因素后，風(fēng)險(xiǎn)評(píng)分是獨(dú)立的預(yù)后指標(biāo)。該標(biāo)志的預(yù)測(cè)價(jià)值在GSE31210數(shù)據(jù)集中得到了驗(yàn)證。為更好地預(yù)測(cè)LUAD患者OS，將風(fēng)險(xiǎn)評(píng)分與臨床特征相結(jié)合，建立了預(yù)測(cè)LUAD患者一、三、五年生存率的列線圖。此外，預(yù)后列線圖校準(zhǔn)曲線顯示，各OS的預(yù)測(cè)存活率和觀察存活率間有較好的一致性。

GO和KEGG富集分析表明，風(fēng)險(xiǎn)評(píng)分改變主要與小分子分解過(guò)程、嘌呤復(fù)合物代謝過(guò)程、線粒體基質(zhì)轉(zhuǎn)移酶復(fù)合物、磷酸酯水解酶活性、輔因子生物合成等通路有關(guān)。已有證據(jù)表明嘌呤類化合物可通過(guò)與katanin相互作用誘導(dǎo)微管斷裂和肺癌細(xì)胞死亡［27］；而另一方面也有研究指出腺嘌呤合成代謝可決定腫瘤干細(xì)胞樣特性和吉非替尼耐藥性產(chǎn)生［28］。輔酶代謝上調(diào)的代謝表達(dá)亞型與腫瘤預(yù)后較差相關(guān)，如被K63-多泛素化后的輔酶DDX17可通過(guò)對(duì)miRNA合成和組蛋白修飾的調(diào)控作用控制腫瘤干細(xì)胞樣特征［29-30］。

基于MRGs的預(yù)后模型可用于評(píng)估潛在的臨床結(jié)局和免疫細(xì)胞浸潤(rùn)。風(fēng)險(xiǎn)評(píng)分和免疫評(píng)分相關(guān)性分析表明，高風(fēng)險(xiǎn)組中NK細(xì)胞、Th1細(xì)胞、Th2細(xì)胞的免疫細(xì)胞亞群顯著上調(diào)；B細(xì)胞、肥大細(xì)胞、中性粒細(xì)胞、T輔助細(xì)胞、iDCs免疫細(xì)胞亞群顯著下調(diào)。Th2細(xì)胞、肥大細(xì)胞和iDCs在GSE84426隊(duì)列中得到了驗(yàn)證。至于免疫相關(guān)通路，APC共刺激、炎癥加劇、Ⅰ型主要組織相容性復(fù)合體和副炎癥顯著上調(diào)，人類白細(xì)胞抗原和Ⅱ型IFN反應(yīng)顯著下調(diào)，副炎癥和Ⅱ型IFN反應(yīng)在GSE84426隊(duì)列中得到了驗(yàn)證。

免疫微環(huán)境會(huì)影響腫瘤發(fā)展、侵襲和轉(zhuǎn)移，不同免疫細(xì)胞浸潤(rùn)程度往往提示不同預(yù)后［31］。已知預(yù)后較差的LUAD患者NK細(xì)胞、Th1細(xì)胞、Th2細(xì)胞、巨噬細(xì)胞、Treg往往顯著上調(diào)，與本研究免疫細(xì)胞浸潤(rùn)程度差異分析結(jié)果相符，且巨噬細(xì)胞、Th2細(xì)胞、Treg上調(diào)在GEO中得到了驗(yàn)證，因此高風(fēng)險(xiǎn)組預(yù)后較低風(fēng)險(xiǎn)組預(yù)后差［20，32-34］。相關(guān)研究表明：活化的NK細(xì)胞可用于抗腫瘤治療，并控制腫瘤生長(zhǎng)和轉(zhuǎn)移擴(kuò)散，但腫瘤微環(huán)境含有的活性氧可抑制NK細(xì)胞抗腫瘤活性，可能導(dǎo)致腫瘤中巨噬細(xì)胞和NK細(xì)胞百分比往往低于正常肺組織［32，35-36］。但LUAD高風(fēng)險(xiǎn)組中，巨噬細(xì)胞和NK細(xì)胞更高的浸潤(rùn)程度卻提示更差的預(yù)后。巨噬細(xì)胞根據(jù)功能活性可分為M1型和M2型，不同LUAD中兩種亞型巨噬細(xì)胞數(shù)量與分布差異顯著，導(dǎo)致LUAD患者不同預(yù)后，研究表明，通常位于腫瘤細(xì)胞池的M1巨噬細(xì)胞與較好的預(yù)后相關(guān)，而在腫瘤間質(zhì)中更為豐富的M2巨噬細(xì)胞與較差預(yù)后相關(guān)［20］。Treg通過(guò)上調(diào)抑制腫瘤相關(guān)第三級(jí)淋巴樣結(jié)構(gòu)中的抗腫瘤反應(yīng)，導(dǎo)致更差預(yù)后［37］。肺癌復(fù)發(fā)患者中，Th2浸潤(rùn)程度顯著上升，Th1/Th2降低。高水平的2型細(xì)胞因子（如IL-10）有助于腫瘤細(xì)胞逃脫免疫監(jiān)視，或許解釋了Th2浸潤(rùn)水平上調(diào)LUAD患者往往預(yù)后不良［37］。

免疫相關(guān)途徑方面，APC共刺激、炎癥加劇、Ⅰ型主要組織相容性復(fù)合物和副炎癥途徑顯著上調(diào)。副炎癥上調(diào)和Ⅱ型IFN反應(yīng)下調(diào)可能與LUAD更好預(yù)后有關(guān)，在TCGA與GEO兩個(gè)隊(duì)列中均得到了證實(shí)：高風(fēng)險(xiǎn)組患者副炎癥途徑上調(diào)而Ⅱ型IFN途徑下調(diào)。程序性定向死亡-1（PD-1-directed）免疫檢查點(diǎn)阻斷可在多種晚期惡性腫瘤中產(chǎn)生持久抗腫瘤活性［38］；Ⅱ型IFN是癌癥和宿主細(xì)胞中PD-L1表達(dá)的關(guān)鍵驅(qū)動(dòng)因素［39］；Ⅱ型IFN反應(yīng)發(fā)生有利于抗腫瘤活性增強(qiáng)，Ⅱ型IFN減少將抑制腫瘤細(xì)胞程序性死亡，因此Ⅱ型IFN途徑下調(diào)預(yù)示不良預(yù)后。

此外，12個(gè)基因（TK1、GCLC、PFKP、CA9、RRM2、HAL、CEL、CA4、CYP2D6、GSTA3、FMO2和ENO3）被確定為關(guān)鍵預(yù)后基因，并可根據(jù)單因素COX分析HR和Lasso模型Coef值確定TK1、GCLC、PFKP、CA9、RRM2和HAL這6個(gè)基因?yàn)槲ｋU(xiǎn)基因，均在高風(fēng)險(xiǎn)LUAD樣本中高表達(dá)，與風(fēng)險(xiǎn)評(píng)分和不良預(yù)后呈正相關(guān)。而CEL、CA4、CYP2D6、GSTA3、FMO2和ENO3這6個(gè)基因被認(rèn)定為保護(hù)基因，其高表達(dá)可降低風(fēng)險(xiǎn)評(píng)分，與LUAD患者不良預(yù)后呈負(fù)相關(guān)。RRM2是核糖核苷酸還原酶催化亞基，可調(diào)節(jié)酶活性，在DNA復(fù)制和修復(fù)中起重要作用，并與細(xì)胞周期、p53信號(hào)傳導(dǎo)途徑、細(xì)胞凋亡等通路相關(guān)。RRM2可獨(dú)立預(yù)測(cè)LUAD預(yù)后，并與免疫浸潤(rùn)有關(guān)，與B細(xì)胞標(biāo)志基因的緊密關(guān)系是免疫應(yīng)答的潛在中心［40］。相關(guān)研究表明：RRM2過(guò)表達(dá)促進(jìn)了Bcl-2和E-鈣黏蛋白信號(hào)通路激活，并抑制p53信號(hào)通路激活，促進(jìn)LUAD細(xì)胞增殖和侵襲［41］。多種途徑（包括糖代謝途徑）可能參與PFKP調(diào)控。有研究從抑制PFKP表達(dá)角度證明：PFKP表達(dá)降低的肺癌細(xì)胞系中葡萄糖攝取率、乳酸水平和三磷酸腺苷濃度顯著降低，可能導(dǎo)致增殖速率、集落形成能力顯著降低和G2/M周期細(xì)胞百分比提高［42］。表明PFKP可能是LUAD的促癌基因。CA4編碼的蛋白是12種人活性同工酶之一，表達(dá)于毛細(xì)血管床、腎臟壁膜、肺微血管和脈絡(luò)膜毛細(xì)血管。已有研究證明CA4低表達(dá)可促進(jìn)癌細(xì)胞增殖，與本研究結(jié)論相同［43］。通過(guò)慢病毒介導(dǎo)技術(shù)建立CA4過(guò)表達(dá)細(xì)胞系證明：CA4過(guò)表達(dá)后，細(xì)胞增殖能力顯著下降，由此得出結(jié)論，CA4可能是抑癌基因，可能與CA4抑制G1期細(xì)胞增殖，誘導(dǎo)細(xì)胞凋亡有關(guān)［44］。與同為碳酸酐酶家族的CA4不同，較少有實(shí)驗(yàn)闡述CA9在LUAD中的致癌機(jī)制，一項(xiàng)對(duì)惡性間皮瘤研究發(fā)現(xiàn)，CA9被特異性抑制劑抑制或在低氧條件下被敲低會(huì)阻礙腫瘤細(xì)胞增殖和遷移，最終導(dǎo)致腫瘤細(xì)胞死亡，證明CA9與腫瘤不良預(yù)后相關(guān)，與本研究結(jié)論一致［45］。表明CA9可能是包括肺癌、惡性間皮瘤等在內(nèi)的多種腫瘤致癌基因。

本研究也存在一定局限性。首先，本研究從公共數(shù)據(jù)庫(kù)中獲得了LUAD患者相應(yīng)資料，無(wú)法排除轉(zhuǎn)錄本可能存在的選擇偏倚。其次，對(duì)模型基因功能的驗(yàn)證需要更多實(shí)驗(yàn)支持。另外，預(yù)測(cè)模型的穩(wěn)定性需在大型前瞻性研究中得到驗(yàn)證。

總之，本研究根據(jù)TCGA和GEO數(shù)據(jù)庫(kù)篩選出的12個(gè)MRGs是LUAD患者生存率、TIME的潛在預(yù)測(cè)指標(biāo)。基于12個(gè)MRGs構(gòu)建的預(yù)測(cè)模型是與OS相關(guān)的獨(dú)立預(yù)測(cè)因素，模型評(píng)分越高，預(yù)后越差?；谠撃Ｐ?，進(jìn)一步開(kāi)發(fā)了一個(gè)方便臨床醫(yī)師預(yù)測(cè)LUAD患者一、三、五年生存率的列線圖，具有良好預(yù)測(cè)性能。此外，本研究還為L(zhǎng)UAD腫瘤微環(huán)境中免疫細(xì)胞浸潤(rùn)研究提供了新的見(jiàn)解，有助于實(shí)現(xiàn)精準(zhǔn)診斷與治療。