基于生物信息學分析胃癌鐵死亡相關(guān)基因預后及免疫細胞浸潤的臨床和功能意義①

羅連響 陳心銘 鄭羽詩 羅 輝 （廣東醫(yī)科大學海洋醫(yī)藥研究院，湛江 524023）

鐵死亡（ferroptosis）是一種區(qū)別于細胞凋亡、細 胞壞死、細胞自噬和細胞焦亡的新型細胞程序性死亡方式，由脂質(zhì)過氧化和致命活性氧積累導致，與腫瘤病理生理過程密切相關(guān)［1-3］。鐵死亡可能在癌癥治療中，尤其是在根除對傳統(tǒng)治療有抵抗力的惡性腫瘤時引發(fā)，癌細胞比正常細胞對鐵的需求量大，這種鐵依賴性可使癌細胞更易受鐵死亡影響。全新的鐵死亡誘導劑有望成為治療耐藥性癌癥的新方法［4-5］。

胃癌（gastric carcinoma，GC）是全球第五大常見癌癥和第三大癌癥死亡原因，是一種分子和表型高度異質(zhì)性的消化道疾病，世界范圍內(nèi)存活率僅為30%，由于人口老齡化，每年新發(fā)病例仍在增加［6-8］。胃癌發(fā)生與幽門螺桿菌感染和年齡等因素相關(guān)［9-10］。隨著放療和化療技術(shù)進步，早期GC 的五年生存率可達到95%，但多數(shù)患者診斷時已經(jīng)是晚期，晚期GC 預后較差，主要治療方法為分子靶向治療、免疫治療及新輔助化療［11-12］。

GC中一些基因在調(diào)節(jié)鐵死亡中起重要作用，如GPX4 和SLC7A11 表達可增加ROS 積累從而促進GC 細胞鐵死亡，SCD1、EIF4A1 和perilipin2 表達可抑制GC 細胞鐵死亡［13-16］。此外，TIMP1 與鐵死亡存在聯(lián)系，在GC中過度表達［17］；TIMP1可能是GC預后的生物標志物［18］；CDO1 對胃癌細胞中Erastin 誘導的鐵死亡具有調(diào)控作用［19］。但鐵死亡相關(guān)基因（fer?roptosis-related gene，F(xiàn)RG）是否與GC 患者預后相關(guān)仍有待研究。

因此，本研究基于TCGA 患者mRNA 表達譜和臨床數(shù)據(jù)構(gòu)建了一個與鐵死亡相關(guān)的預后模型，在GEO 隊列中進行了驗證，并對其功能富集和免疫打分相關(guān)性進行分析，旨在開發(fā)一種預測GC 患者總生存率（overall survival，OS）的FRG簽名和列線圖模型，探討其可能作用機制，為GC 患者個性化治療提供參考。

1 材料與方法

1.1 表達數(shù)據(jù)和臨床數(shù)據(jù)收集 研究流程圖如圖1所示。通過TCGA（The Cancer Genome Atlas）數(shù)據(jù)庫獲取胃癌（Stomach adenocarcinoma，STAD）樣本基因表達及相關(guān)臨床信息，共得到403例STAD樣本用于訓練組研究。通過美國國立生物信息中心（National Center for Biotechnology Information，NCBI）的GEO（Gene Expression Omnibus）數(shù)據(jù)庫獲取數(shù)據(jù)集GSE84426的基因表達數(shù)據(jù)和臨床信息，共得到76例GC 樣本被納入驗證組研究。兩組臨床特征如表1所示。胃腺癌代表了多數(shù)胃惡性腫瘤［20］，占胃惡性腫瘤的95%，因此來自TCGA-STAD 和GSE84426 的樣本可代表胃腺癌樣本。

圖1 FRG預后模型構(gòu)建Fig.1 Construction of prognosis model of FRG

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

1.2 篩選預后相關(guān)的差異基因 從FerrDb（http：//www.zhounan.org/ferrdb/）收集了290個FRG。采用R軟件limma 包對數(shù)據(jù)進行整理，提取出245FRG，篩選出115 個差異表達基因（differential expression genes，DEG,fdr<0.05，|logFC|>0.5）。采用R 軟件survival 包、單因素Cox 分析得到28個預后相關(guān)基因（P<0.05），比較后得到的10個差異表達的FRG具有預后價值用于構(gòu)建預后模型。

1.3 預后模型的構(gòu)建和驗證 基于上述10 個基因，采用R 軟件glmnet包進行LASSO-Cox回歸分析，構(gòu)建預后模型?；颊唢L險評分=∑（每個基因表達水平×相應系數(shù)）。根據(jù)風險評分的中位數(shù)將患者分為高風險組和低風險組。受試者操作特征（ROC）分析檢驗獨立風險因素預測生存率的敏感性和特異性，評估模型預測準確性。t-SNE 和PCA 分析探討高風險組和低風險組分布。

1.4 功能富集分析 采用R 軟件clusterProfiler 包進行GO（Gene Ontology）和KEGG（Kyoto Encyclopedia of Genes and Genomes）分析，探討不同分子機制及高、低風險患者差異。GO 分析包括生物學過程、細胞成分和分子功能。

1.5 免疫細胞和免疫相關(guān)通路評分的相關(guān)性分析 采用R 軟件GSVA 包中的單樣本基因集富集分析（ssGSEA）進一步量化16 個免疫細胞的浸潤水平和13個免疫相關(guān)通路活性。

1.6 統(tǒng)計學分析t檢驗用于確定腫瘤組織和正常組織中差異表達的FRG。采用卡方檢驗比較各組比例構(gòu)成差異。Mann-Whitney 檢驗評估風險組間免疫細胞或途徑ssGSEA 得分差異，并采用BH 方法調(diào)整P值。采用Kaplan-Meier分析和對數(shù)秩檢驗（Logrank test）比較各組OS。單因素和多因素Cox 回歸分析確定OS 的獨立預后因素。所有統(tǒng)計分析均采用R 軟件（4.0.3）完成。P<0.05 為差異具有統(tǒng)計學意義。

2 結(jié)果

2.1 TCGA-STAD隊列中預后相關(guān)DEG篩選 TCGASTAD 隊列中篩選出115 個鐵死亡相關(guān)DEG，其中14 個與OS有關(guān)（圖2A）。14 個鐵死亡相關(guān)DEG 中，8個基因在腫瘤組織中上調(diào)，6個基因在腫瘤組織中下調(diào)（圖2B）。單因素Cox 分析顯示，14 個基因均與GC患者OS顯著相關(guān)，其中9個為危險基因（HR>1），5 個為保護基因（HR<1）（圖2C）。相關(guān)性網(wǎng)絡顯示14 個有預后價值的鐵死亡相關(guān)DEG 相關(guān)（圖2D）。STRING 數(shù)據(jù)庫下載的PPI 網(wǎng)絡表明基因間相互作用（圖2E）。14 個基因中ATG16L1、CAV1、CDO1、DUSP1、GABARAPL1、HIF1A、LONP1、MYB、NFE2L2、NOX4、SLC1A5和TGFBR1在腫瘤組織和正常組織中表達差異有統(tǒng)計學意義（P<0.05，圖2F），其中ATG16L1、HIF1A、LONP1、MYB、NOX4、SLC1A5 和TGFBR1 在腫瘤組織表達上調(diào)，為致癌基因；CAV1、CDO1、DUSP1、GABARAPL1 和NFE2L2 在腫瘤組織表達下調(diào)，為抑癌基因（圖2F）。

圖2 TCGA-STAD隊列中有預后價值的DEG篩選Fig.2 Screening of DEG with prognostic value in TCGASTAD cohort

2.2 10 基因預后模型構(gòu)建 LASSO-Cox 回歸分析得到10 個有預后價值的DEG（表2，基因全稱來自https：//www. 360zhyx. com/），構(gòu)建預后模型（圖3）。以風險評分中位數(shù)為臨界值，將TCGA-STAD 隊列患者分為高風險組（n=159）和低風險組（n=158）。Kaplan-Meier曲線顯示高風險組生存時間縮短（圖4A，P<0.01）。ROC 曲線評價模型預測性能，曲線下面積（AUC）在3 年時達到0.693，4 年時達到0.695，5 年時達到0.700（圖4B）。圖4C 顯示患者被分為高、低風險兩組。高風險組患者生存時間較低風險患者生存時間縮短（圖4D）。PCA 和t-SNE 分析顯示，不同風險組患者分為兩個方向（圖4E、F）。

表2 預后簽名中的基因Tab.2 Genes in prognostic gene signatures

圖3 TCGA-STAD隊列10基因預后模型構(gòu)建Fig.3 Construction of 10 gene prognosis model in TCGASTAD cohort

圖4 TCGA-STAD隊列基于10基因預后模型的生存分析Fig.4 Survival analysis in TCGA-STAD cohort based on 10 gene prognostic model

2.3 10 基因預后模型的驗證 將GSE84426 隊列患者分為高風險組（n=38）和低風險組（n=38）。GSE84426隊列生存分析結(jié)果與TCGA-STAD隊列生存分析結(jié)果相近（圖5，P<0.05，AUC 在3 年時達到0.678，4年時達到0.686，5年時達到0.691）。

圖5 預后模型驗證Fig.5 Validation of prognostic model

2.4 獨立預后分析 單因素Cox 回歸分析中，分期較高是TCGA-STAD 隊列和GSE84426 隊列患者OS的不利因素（TCGA-STAD：HR=1.833，95%CI=1.256～2.674，P=0.002；GSE84426：HR=2.690，95%CI=1.042～6.943，P=0.041，圖6A、C）。多因素Cox回歸分析中，年齡較大在TCGA-STAD 隊列中顯示出更差的OS（HR=2.038，95%CI=1.395～2.976，P<0.001，圖6B）。此外，高風險是兩組患者OS 的不利因素（TCGA-BRCA：HR=4.238，95%CI=2.507～7.613，P<0.001；GSE84426：HR=2.137，95%CI=1.057～4.319，P=0.034，圖6B、D）。因此，風險評分是OS 的獨立預后因素。圖6E顯示了TCGA-STAD 隊列構(gòu)建的列線圖，采用4個預后因素（年齡、性別、分級和分期），可根據(jù)影響生存的預后因素的比例對患者進行打分。

圖6 獨立預后分析Fig.6 Independent prognostic analysis

2.5 功能富集分析 GO 分析顯示，來自TCGA 隊列風險組間的DEG 主要富集于細胞外基質(zhì)組織和細胞外結(jié)構(gòu)組織。細胞成分方面，DEG 主要富集于含膠原細胞外基質(zhì)（圖7A）。KEGG 途徑分析顯示，黏著斑和PI3K-Akt信號通路顯著富集（圖7B）。

圖7 GO和KEGG功能富集分析Fig.7 Enrichment analysis of GO and KEGG functions

2.6 免疫細胞和免疫相關(guān)途徑 TCGA 隊列高風險組中，B 細胞、CD8+T 細胞、樹突狀細胞、未成熟樹突狀細胞、巨噬細胞、肥大細胞、中性粒細胞、漿細胞樣樹突狀細胞、輔助性T 細胞、濾泡輔助性T 細胞、TIL細胞和調(diào)節(jié)性T細胞的免疫細胞亞群顯著上調(diào)（所有調(diào)整后的P<0.05，圖8A）。至于免疫相關(guān)通路，APC 共刺激、CCR、檢查點、副炎癥、T 細胞共刺激和Ⅱ型IFN 反應顯著上調(diào)（均經(jīng)校正P<0.05，圖8B）。通過GSE84426 隊列的高風險組驗證，B 細胞、濾泡輔助性T 細胞、TIL 細胞、T 細胞共刺激和Ⅱ型IFN反應均被證實（均經(jīng)調(diào)整P<0.05，圖8C、D）。

圖8 TCGA-STAD 隊列和GSE84426 隊列高風險組和低風險組ssGSEA得分Fig.8 ssGSEA scores of high-risk group and low-risk group in TCGA-STAD cohort and GSE84426 cohort

3 討論

GC患者晚期預后較差，五年生存率較低使開發(fā)穩(wěn)定的預后指標變得非常重要。本研究基于TCGA數(shù)據(jù)建立了一個包含10 個基因簽名的鐵死亡相關(guān)模型預測GC 預后，并在GSE84426 隊列中驗證其預測準確性。上述兩個隊列中，高風險組GC 患者生存期均短于低風險組。此外，風險評分是GC 患者獨立的預測因子，與臨床特征和免疫功能密切相關(guān)。

本研究通過對10 個鐵死亡相關(guān)基因的功能分析，確 定MYB、SLC1A5、NFE2L2、ATG16L1 和LONP1 為 保 護 基 因，TGFBR1、ZFP36、DUSP1、HIF1A和NOX4為危險基因。根據(jù)FerrDb，TGFBR1、HIF1A、MYB、SLC1A5、DUSP1 和NOX4 是促進鐵死亡的驅(qū)動因素，而NFE2L2 和ZFP36 是預防鐵死亡的抑制劑（http：//www. zhounan. org/ferrdb/）。對于TGFBR1，Linc00462過表達提高TGFBR1和TGFBR2表達，從而促進胰腺癌侵入性［21-22］。此外，TGFBR1表達可促進GC 發(fā)展［23］。MYB 在多種白血病和癌癥中過度表達，與ROS 遞質(zhì)改變有關(guān)［24］。SLC1A5 表達的蛋白是一種細胞表面轉(zhuǎn)運蛋白，轉(zhuǎn)運中性氨基酸［25］。miR-137 通過直接靶向谷氨酰胺轉(zhuǎn)運體SLC1A5 負性調(diào)節(jié)黑色素瘤細胞鐵死亡［26］；SLC1A5表達可抑制GC 細胞增殖［27］。NFE2L2（NRF2）可促進癌細胞增殖，其表達可能增強GC 細胞對治療的抵抗力［28-29］。另外，NFE2L2 過表達促進了對鐵死亡的抵抗。ARF 介導的鐵死亡在很大程度上被NFE2L2 共表達所消除。NFE2L2 基因抑制增強了鐵死亡易感性（http：//www.zhounan.org/ferrdb/）。研究發(fā)現(xiàn)ZFP36在肝星狀細胞鐵死亡調(diào)控中起關(guān)鍵作用，可抑制鐵死亡，并確定ZFP36自噬依賴性鐵死亡是治療肝纖維化的一個潛在靶點，ZFP36 質(zhì)粒通過破壞ATG16L1 mRNA 穩(wěn)定性抑制自噬激活，ATG16L1 質(zhì)粒可消除ZFP36 質(zhì)粒對鐵死亡的抑制作用，與GC 易感性相關(guān)［30-31］。ATG16L1是自噬相關(guān)的基因［32］，本研究發(fā)現(xiàn)ATG16L1 與鐵死亡有關(guān)，其在GC 組織中表達上調(diào)，且與正常組織表達差異顯著。有研究闡明了LONP1 與鐵死亡的關(guān)系，LONP1抑制其對胰腺導管腺癌PANC1 細胞鐵死亡的保護作用［33］。DUSP1 在erastin 或RSL3 誘導的鐵死亡過程中表達上調(diào)（http：//www.zhounan.org/ferrdb/）不穩(wěn)定的HIF1A 促進腫瘤細胞鐵死亡［34］。土荊皮乙酸通過激活NOX4 和抑制xCT 可引起膠質(zhì)瘤細胞鐵死亡［35］。

單因素Cox 回歸分析中，分期和風險評分與OS顯著相關(guān)。多因素Cox 回歸分析顯示，風險評分是一個獨立的預后因素。為更好地預測GC 患者的OS，本研究將風險評分與臨床特征相結(jié)合，建立了GC患者列線圖。

功能富集分析表明，風險評分改變主要與細胞外基質(zhì)組織、細胞外結(jié)構(gòu)組織和PI3K-Akt 信號通路有關(guān)。PI3K-Akt 信號通路參與細胞周期進展、凋亡和腫瘤轉(zhuǎn)化，是腫瘤治療靶點，與較差預后有關(guān)。作為正常細胞過程中最重要的信號通路之一，PI3K-Akt 信號通路可作為治療GC 和胰腺癌的潛在治療策略［36-37］。另外，PI3K 的活動也與其他腫瘤相關(guān)，包括肺癌、黑色素瘤和白血?。?8］。

免疫評分和免疫評分相關(guān)性分析表明，高風險組中，B 細胞、CD8+T 細胞、樹突狀細胞、未成熟樹突狀細胞、巨噬細胞、肥大細胞、中性粒細胞、漿細胞樣樹突狀細胞、輔助性T 細胞、濾泡輔助性T 細胞、TIL細胞和調(diào)節(jié)性T細胞的免疫細胞亞群顯著上調(diào)，B 細胞、濾泡輔助性T 細胞、TIL 細胞在GSE84426 隊列中得到了驗證。已知B 細胞、T 細胞、巨噬細胞、肥大細胞、中性粒細胞、TIL 細胞和調(diào)節(jié)性T 細胞的增加預示GC 患者不良預后，因此高風險組預后較低風險組預后差［39-44］。在免疫相關(guān)途徑方面，APC共刺激、CCR、檢查點、副炎癥、T 細胞共刺激和Ⅱ型IFN 反應顯著上調(diào)，T 細胞共刺激和Ⅱ型IFN 反應被驗證組證實，均為上調(diào)。程序性死亡-1-導向（PD-1-導向）免疫檢查點阻斷法可將持久的抗腫瘤活性賦予多種晚期惡性腫瘤。Ⅱ型IFN 反應是癌細胞和宿主細胞中PD-L1表達的關(guān)鍵驅(qū)動因素［45］。高風險組Ⅱ型IFN 反應在這兩個隊列中上調(diào)，而Ⅱ型IFN 反應有利于抗腫瘤活性增強，因此高風險組預后較差的原因可能是T 細胞共刺激的作用大于Ⅱ型IFN 反應的作用。本研究也有局限性，如僅采用公共數(shù)據(jù)庫，需要更多臨床驗證。

總之，本研究應用TCGA 和GEO 數(shù)據(jù)集構(gòu)建并驗證了鐵死亡相關(guān)基因的預后模型，該模型是與OS相關(guān)的獨立因素。此外，本研究還開發(fā)了一個列線圖，醫(yī)師可應用列線圖提高高風險GC 患者識別的準確性，實現(xiàn)精準治療。