王楊淑怡 葉 威 林志豪 黃約諾 方 濤 董曉亭

（浙江中醫(yī)藥大學(xué)附屬溫州中醫(yī)院，浙江 溫州 325000）

細(xì)胞焦亡是一種程序性細(xì)胞死亡的形式，由炎癥小體觸發(fā)，被半胱氨酸天冬氨酸特異性蛋白酶（cysteine-containing aspartate-specific protease，Caspase）執(zhí)行，以質(zhì)膜的快速破裂和促炎因子的釋放為特征[1]。細(xì)胞焦亡可分為Caspase-1依賴的經(jīng)典途徑和Caspase-4、Caspase-5、Caspase-11依賴的非經(jīng)典途徑[1]。Caspase-1依賴的經(jīng)典途徑可以被細(xì)菌、病毒等[2]激活，使消皮素D（gasdermin D，GSDMD）家族蛋白被剪切生成N端和C端，GSDMD與細(xì)胞膜磷脂結(jié)合形成10～20 nm的小孔，細(xì)胞內(nèi)容物通過孔隙緩慢釋放，觸發(fā)炎癥反應(yīng)[3]。細(xì)胞逐漸變平，產(chǎn)生 1～5 μm的凋亡囊泡狀突起，并逐漸膨脹，直至質(zhì)膜破裂[3]。Caspase-4、Caspase-5、Caspase-11依賴的非經(jīng)典途徑可由脂多糖等致炎因子激活，生成GSDMD氮端活性域的肽段，引起炎癥反應(yīng)，啟動細(xì)胞焦亡[3]。

越來越多的研究發(fā)現(xiàn)，細(xì)胞焦亡與腫瘤的發(fā)生、發(fā)展密切相關(guān)[4-5]。長期、慢性的炎癥可能導(dǎo)致局部組織的異常增生，進(jìn)而發(fā)生癌變[4-5]。GAO等[2]的研究結(jié)果顯示，非小細(xì)胞肺癌（non-small cell lung cancer，NSCLC）患者GSDMD表達(dá)顯著上調(diào)，參與了腫瘤的病理進(jìn)展。WANG等[6]發(fā)現(xiàn)，辛伐他丁通過激活炎癥小體核苷酸結(jié)合寡聚化結(jié)構(gòu)域樣受體3（nucleotidebinding oligomerization domain-like receptor 3，NLRP3）激活Caspase-1，介導(dǎo)細(xì)胞焦亡的經(jīng)典途徑，抑制癌細(xì)胞的增殖和遷移。肺腺癌（lung adenocarcinoma，LUAD）是NSCLC的常見病理類型，患者的生存率和預(yù)后往往較差，因此制定有效的LUAD治療策略迫在眉睫[7]。目前，LUAD與細(xì)胞焦亡之間關(guān)系的研究較少。本研究通過機(jī)器學(xué)習(xí)算法分析細(xì)胞焦亡基因在LUAD癌組織中的表達(dá)特征，及其在患者預(yù)后評估中的價值，為制定LUAD的干預(yù)措施提供一個新的思路。

1 材料和方法

1.1 數(shù)據(jù)收集

從腫瘤基因組圖譜（the Cancer Genome Atlas，TCGA）數(shù)據(jù)庫（https：//portal.gdc.cancer.gov/）下載LUAD患者的轉(zhuǎn)錄組數(shù)據(jù)和對應(yīng)的臨床資料（截至2021月1月28日），包括535例LUAD癌組織樣本和59例癌旁組織樣本。TCGA-LUAD數(shù)據(jù)集采用R軟件limma程序包進(jìn)行標(biāo)準(zhǔn)化處理（TCGA隊列）。

從基因表達(dá)綜合（Gene Expression Omnibus，GEO）數(shù)據(jù)庫（http：//www.ncbi.nlm.nih.gov/geo）下載符合條件的LUAD轉(zhuǎn)錄組數(shù)據(jù)和對應(yīng)的臨床數(shù)據(jù)。篩選條件：1）樣本量較大；2）包含年齡、性別、生存狀態(tài)、生存期等臨床信息。根據(jù)篩選條件納入符合條件的GEO隊列（GSE72094數(shù)據(jù)集），并進(jìn)行注釋。

1.2 差異表達(dá)基因分析

從GeneGards數(shù)據(jù)庫中鑒定出r>0.1的細(xì)胞焦亡相關(guān)基因。在Pubmed數(shù)據(jù)庫（https：//pubmed.ncbi.nlm.nih.gov/）中檢索關(guān)鍵詞“pyroptosis”，驗(yàn)證上述焦亡相關(guān)基因，共鑒定出52個細(xì)胞焦亡基因。通過limma程序包提取TCGA焦亡基因的表達(dá)量，將差異倍數(shù)（fold change，F(xiàn)C）≥1，即|log2FC|≥1且P<0.05設(shè)置為過濾條件，篩選出聚類分型間的差異表達(dá)基因。鑒定出LUAD癌組織和癌旁組織間的差異表達(dá)基因，采用pheatmap程序包進(jìn)行可視化。通過STRING數(shù)據(jù)庫（https：//www.string-db.org/）獲取差異表達(dá)基因的蛋白質(zhì)-蛋白質(zhì)互作（proteinprotein interaction，PPI）網(wǎng)絡(luò)。運(yùn)用corrplot程序包計算差異表達(dá)基因的r值，以r=0.5為閾值進(jìn)行篩選。

1.3 差異表達(dá)基因的共識聚類分析

將樣本的生存數(shù)據(jù)和樣本（含有差異表達(dá)焦亡基因的表達(dá)量）進(jìn)行合并、標(biāo)準(zhǔn)化等處理。利用survival程序包對基因進(jìn)行Cox回歸分析，篩選與生存相關(guān)的焦亡基因。根據(jù)生存相關(guān)焦亡基因的表達(dá)特征，采用limma程序包和ConsensusClusterPlus程序包對樣本進(jìn)行聚類分析，分別獲取條件滿足組間關(guān)系不緊密和組間關(guān)系緊密的焦亡基因聚類分型。采用survival程序包和survminer程序包對聚類分型繪制Kaplan-Meier生存曲線，用秩和法進(jìn)行檢驗(yàn)。采用pheatmap程序包繪制熱圖，展現(xiàn)不同聚類分型在臨床特征中的分布情況。

1.4 構(gòu)建基于LUAD細(xì)胞焦亡基因的預(yù)后風(fēng)險模型

利用limma程序包對聚類分型進(jìn)行差異分析，將|log2FC|≥1且P<0.05設(shè)置為過濾條件，篩選出焦亡模型間的差異表達(dá)基因。采用Cox回歸分析評估每個差異表達(dá)基因與LUAD患者生存的關(guān)系。以P<0.05為篩選標(biāo)準(zhǔn)，鑒定出生存相關(guān)差異表達(dá)基因。采用glmnet程序包進(jìn)行最小絕對收縮和選擇算子（least absolute shrinkage and selection operator，LASSO）Cox回歸分析，縮小候選基因的范圍，建立預(yù)后風(fēng)險模型。懲罰系數(shù)（λ）是由最小參數(shù)決定候選基因及其系數(shù)。最后基于所有候選基因的系數(shù)值與基因表達(dá)總和的積確定風(fēng)險評分。

將TCGA隊列作為測試集，用于內(nèi)部評估；GEO隊列作為驗(yàn)證集，用于外部評估。根據(jù)風(fēng)險評分的中位值將樣本劃分為高風(fēng)險組、低風(fēng)險組。采用timeROC程序包進(jìn)行1年生存、3年生存和5年生存的受試者工作特征（receiver operating characteristic，ROC）曲線分析。此外，根據(jù)TCGA風(fēng)險評分中位值將GEO隊列分為高風(fēng)險組、低風(fēng)險組，采用timeROC包繪制1年生存、3年生存和5年生存的ROC曲線進(jìn)行驗(yàn)證，曲線下面積（area under curve，AUC）>0.5提示模型有較好的預(yù)測價值。

1.5 風(fēng)險評分和LUAD患者生存之間的關(guān)系

采用χ2檢驗(yàn)分析測試集中高風(fēng)險組、低風(fēng)險組與LUAD患者臨床病理特征之間的相關(guān)性，以熱圖形式呈現(xiàn)。采用Kaplan-Meier生存曲線分別比較測試集和驗(yàn)證集不同風(fēng)險組之間生存情況的差異。采用單變量和多變量Cox回歸分析確定測試集和驗(yàn)證集風(fēng)險評分的獨(dú)立性。上述分析過程均采用survival程序包和survminer程序包完成。

1.6 生物學(xué)功能分析

采用limma程序包比較TCGA隊列高風(fēng)險組和低風(fēng)險組之間的差異表達(dá)基因，以|log2FC|≥1且P<0.05為篩選條件。采用clusterProfiler程序包進(jìn)行基因本體論（Gene Ontology，GO）和京都基因與基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）分析。以FDR<0.05為顯著富集的生物學(xué)功能和信號通路。

1.7 免疫浸潤分析

采用GSVA程序包和GSEABase程序包對TCGA隊列和GEO隊列進(jìn)行單樣本基因集富集分析（single sample gene set enrichment analysis，ssGSEA），量化免疫細(xì)胞和免疫相關(guān)通路，采用Mann-Whitney檢驗(yàn)比較TCGA隊列免疫浸潤細(xì)胞和免疫相關(guān)通路的差異。采用CIBERSORT程序包（基于線性支持向量回歸的原理）對人類22種免疫細(xì)胞亞型和免疫功能進(jìn)行量化分析。采用Wilcoxon秩和檢驗(yàn)比較高風(fēng)險組和低風(fēng)險組之間免疫浸潤細(xì)胞的差異。以P<0.05為差異有統(tǒng)計學(xué)意義。

1.8 統(tǒng)計學(xué)分析

采用R 4.0.3軟件進(jìn)行統(tǒng)計分析。

2 結(jié)果

2.1 差異表達(dá)的焦亡基因

共篩選出41個差異表達(dá)的焦亡基因，其中表達(dá)下調(diào)10個（IL6、NLRC4、CASP5、IL1B、NLRP3、IL1A、CASP1、IRF1、CHMP3、NLRP1）、表達(dá)上調(diào)31個（TNF、SCAF11、CHMP2B、ELANE、CHMP6、CASP9、NLRP6、CHMP7、TIRAP、CHMP2A、PLCG1、GSDMD、CASP4、BAX、GPX4、CASP8、CHMP4A、CHMP4B、GSDME、TP53、PJVK、CYCS、CASP3、BAK1、CASP6、CHMP4C、GSDMA、GSDMB、NLRP7、GSDMC、AIM2）（P<0.05），見圖1（a）。PPI分析結(jié)果表明（相互作用分?jǐn)?shù)為0.9），CHMP3、CHMP4C、CHMP2A、CHMP7、CHMP4B、CHMP2B、CHMP6、CHMP4A、CASP6、CASP3、CASP8是核心基因，見圖1（b）。以閾值r=0.5為臨界點(diǎn)繪制相關(guān)網(wǎng)絡(luò)圖，見圖1（c）。

2.2 基于差異表達(dá)基因的腫瘤分型

根據(jù)焦亡基因的表達(dá)特征將LUAD樣本通過聚類分析分型。聚類的變量用字母k表示，同時根據(jù)繪制圖形獲得最合適的分型。結(jié)果顯示，在k=2的情況下，LUAD患者可以分為2個完全不同的、不重疊的聚類分型，見圖2（a）?；虮磉_(dá)譜和臨床病理特征（包括性別、年齡、Stage分期和TMN分期）在聚類分型中的分布見圖2（b）。分型2的生存期顯著長于分型1（P<0.05），見圖2（c）。

圖2 基于差異表達(dá)基因特征的腫瘤分型

2.3 基于LUAD焦亡基因的預(yù)后模型構(gòu)建及其臨床價值

篩選出聚類分型間的差異表達(dá)基因，采用Cox回歸分析鑒定出11個與生存相關(guān)的差異表達(dá)基因，分別為：ANO1、PKHD1L1、FMO2、TDRD1、TBX4、ABCC12、AQP4、CPXM2、WFIKKN2、ATP1A4、IGSF10。風(fēng)險比（hazard ratio，HR）>1的基因?yàn)轱L(fēng)險基因，即基因的表達(dá)量越多，預(yù)后越差，HR<1的基因?yàn)楸Ｗo(hù)基因，即基因的表達(dá)量越多，預(yù)后越好。采用LASSO-Cox回歸分析，并根據(jù)最佳λ值構(gòu)建預(yù)后風(fēng)險模型。Kaplan-Meier生存曲線分析結(jié)果顯示，測試集（TCGA隊列）高風(fēng)險組的生存期顯著短于低風(fēng)險組（P<0.05）。與低風(fēng)險組相比，高風(fēng)險組的死亡例數(shù)更多、生存期更短。ROC曲線分析結(jié)果顯示，預(yù)后風(fēng)險模型判斷TCGA隊列LUAD患者1年生存的AUC為0.744，3年生存的AUC為0.697，5年生存的AUC為0.667。見圖3。

圖3 預(yù)后模型對TCGA隊列的預(yù)后評估效能

Kaplan-Meier生存曲線分析結(jié)果顯示，驗(yàn)證集（CEO隊列）高風(fēng)險組生存期短于低風(fēng)險組（P<0.05）。與低風(fēng)險組比較，高風(fēng)險組生存期更短，死亡率更高。ROC曲線分析結(jié)果顯示，預(yù)后模型判斷GEO隊列LUAD患者1年生存的AUC為0.561、3年生存的AUC為0.611、5年生存的AUC為0.734。見圖4。

圖4 預(yù)后模型對GEO隊列的預(yù)后評估效能

2.4 預(yù)后風(fēng)險模型與LUAD患者預(yù)后的關(guān)系

單因素Cox回歸分析結(jié)果顯示，依據(jù)預(yù)后風(fēng)險模型得出的風(fēng)險評分是TCGA隊列LUAD患者預(yù)后的危險因素[HR=2.911，95%可信區(qū)間（confidence interval，CI）為1.355～6.256，P=0.006]。多因素Cox回歸分析結(jié)果顯示，在排除其他混雜因素（年齡、性別、T分期、N分期、M分期）后，依據(jù)預(yù)后風(fēng)險模型得出的風(fēng)險評分是TCGA隊列LUAD患者預(yù)后的獨(dú)立危險因素（HR=4.815，95%CI為2.963～7.824，P<0.001）。但基于GEO隊列得出了不同的結(jié)果（HR=1.857，95%CI為0.857～4.071，P=0.122）?；赥CGA隊列11個差異表達(dá)基因的熱圖分析結(jié)果顯示，高風(fēng)險組、低風(fēng)險組之間性別、臨床分期、T分期、N分期差異有統(tǒng)計學(xué)意義（P<0.05）。見圖5。

圖5 預(yù)后風(fēng)險模型的Cox回歸分析和TCGA隊列高風(fēng)險組、低風(fēng)險組臨床病理特征的差異熱圖

2.5 TCGA隊列高風(fēng)險組和低風(fēng)險組差異表達(dá)基因的生物學(xué)功能分析

KEGG富集分析和GO富集分析結(jié)果顯示，TCGA隊列高風(fēng)險組和低風(fēng)險組差異表達(dá)基因主要與趨化因子介導(dǎo)腫瘤相關(guān)信號通路、免疫反應(yīng)、細(xì)胞膜功能相關(guān)。見圖6。

圖6 TCGA隊列高風(fēng)險組、低風(fēng)險組差異表達(dá)基因的功能分析和ssGSEA結(jié)果

采用ssGSEA量化TCGA隊列LUAD患者的免疫相關(guān)通路。與高風(fēng)險組比較，低風(fēng)險組B淋巴細(xì)胞、抗原遞呈細(xì)胞、T淋巴細(xì)胞、人類白細(xì)胞抗原、中性粒細(xì)胞、腫瘤浸潤淋巴細(xì)胞、調(diào)節(jié)性T細(xì)胞的比例更高。見圖6。

采用CIBERSORT量化免疫浸潤細(xì)胞。與高風(fēng)險組比較，低風(fēng)險組CD4+記憶T淋巴細(xì)胞、肥大細(xì)胞、樹突狀細(xì)胞比例更高（P<0.05）。見圖6。

3 討論

細(xì)胞焦亡是一種新的促炎程序性細(xì)胞死亡方式，在腫瘤的發(fā)生、發(fā)展中發(fā)揮著雙重作用[1，3]。一方面，在細(xì)胞焦亡的進(jìn)程中釋放大量的炎癥因子，促進(jìn)腫瘤的發(fā)生、發(fā)展[8-9]；另一方面，細(xì)胞焦亡的促炎作用可減輕耐藥性，更好地發(fā)揮藥物的抗癌作用[8-9]。細(xì)胞焦亡可能是一種潛在的腫瘤治療新策略。本研究發(fā)現(xiàn)了41個在LUAD患者癌組織和癌旁組織中呈差異表達(dá)的焦亡基因，其中CHMP3、CHMP4C、CHMP2A、CHMP7、CHMP4B、CHMP2B、CHMP6、CHMP4A、CASP6、CASP3、CASP8在LUAD的發(fā)生、發(fā)展中起重要作用。基于差異基因的表達(dá)特征，本研究將LUAD樣本分成2個聚類分型，與LUAD患者預(yù)后相關(guān)。進(jìn)一步篩選出聚類分型間11個生存相關(guān)的差異表達(dá)基因，包括ANO1、PKHD1L1、FMO2、TDRD1、TBX4、ABCC12、AQP4、CPXM2、WFIKKN2、ATP1A4、IGSF10。生存分析結(jié)果顯示，高表達(dá)的ANO1與較差的預(yù)后相關(guān)；PKHD1L1、FMO2、TDRD1、TBX4、ABCC12、AQP4、CPXM2、WFIKKN2、ATP1A4、IGSF10為保護(hù)基因，其高表達(dá)與較好的預(yù)后相關(guān)。這些基因可能是LUAD和細(xì)胞焦亡聯(lián)系的關(guān)鍵節(jié)點(diǎn)，在疾病進(jìn)展中發(fā)揮著截然不同的作用。通過LASSOCox回歸分析建立了一個焦亡相關(guān)的預(yù)后模型，具有不同的性別、臨床分期、T分期、N分期和M分期的特征。外部驗(yàn)證集獲得相似的結(jié)果，并發(fā)現(xiàn)高、低風(fēng)險組間的差異基因與腫瘤免疫相關(guān)。

本研究KEGG富集分析結(jié)果顯示，差異表達(dá)基因在腫瘤相關(guān)信號通路上存在顯著差異。有研究證實(shí)，細(xì)胞焦亡通過調(diào)控某個靶點(diǎn)或某些信號通路的活性來發(fā)揮抗惡性腫瘤的作用[10]。核因子-κB（nuclear factor-kappa B，NF-κB）信號通路的持續(xù)異常激活能導(dǎo)致腫瘤細(xì)胞在抗凋亡的同時促使其浸潤、轉(zhuǎn)移，在惡性腫瘤的發(fā)生、發(fā)展中起重要作用，抑制NF-κB信號通路的異常激活能夠促進(jìn)包括細(xì)胞焦亡在內(nèi)的程序性死亡，進(jìn)而發(fā)揮抗腫瘤作用[11]。CIBERSORT程序包分析結(jié)果顯示，TCGA隊列高風(fēng)險組與低風(fēng)險組之間免疫浸潤細(xì)胞比例存在差異。細(xì)胞焦亡能調(diào)節(jié)免疫浸潤細(xì)胞的狀態(tài)，調(diào)控腫瘤微環(huán)境的組分。部分腫瘤細(xì)胞發(fā)生焦亡后，釋放的炎性介質(zhì)可激活T細(xì)胞，調(diào)節(jié)腫瘤免疫微環(huán)境，激活機(jī)體的免疫反應(yīng)，發(fā)揮抗腫瘤的功能[12]。

本研究采用ssGSEA、CIBERSORT分析高風(fēng)險組、低風(fēng)險組之間的免疫特性，結(jié)果顯示，低風(fēng)險組擁有更高的T淋巴細(xì)胞、樹突狀細(xì)胞、肥大細(xì)胞，與免疫相關(guān)通路有關(guān)，有著更好的預(yù)后。細(xì)胞焦亡在腫瘤微環(huán)境中能夠激活抗腫瘤免疫，抑制腫瘤生長。HSU等[13]的研究結(jié)果顯示，在GSDME-鼠的三陰性乳腺癌細(xì)胞系EMT6和結(jié)直腸癌細(xì)胞系CT26的微環(huán)境中，CD8+T淋巴細(xì)胞和自然殺傷細(xì)胞浸潤的減少伴隨著腫瘤浸潤淋巴細(xì)胞釋放顆粒酶B和穿孔素的減少[13]。腫瘤細(xì)胞發(fā)生細(xì)胞焦亡時會產(chǎn)生促炎因子和免疫抗原，激活Gasdermin家族蛋白，N端分子結(jié)構(gòu)域在細(xì)胞膜表面聚集形成寡肽空，促進(jìn)細(xì)胞裂解，誘導(dǎo)機(jī)體發(fā)生瀑布級聯(lián)炎癥反應(yīng)，激活以T細(xì)胞為主的抗腫瘤免疫反應(yīng)[12]。

近年來，相關(guān)研究構(gòu)建了很多LUAD的風(fēng)險預(yù)后模型，為LUAD的診斷和預(yù)后提供了新的思路。本研究證實(shí)了焦亡基因?qū)UAD的診斷和預(yù)后具有重要價值，通過LASSO-Cox回歸分析建立風(fēng)險預(yù)后模型，采用傳統(tǒng)方法證明風(fēng)險預(yù)后模型的有效性。細(xì)胞焦亡會激活機(jī)體免疫反應(yīng)，發(fā)揮抗癌的作用，設(shè)計誘導(dǎo)細(xì)胞焦亡的治療策略可以加強(qiáng)抗腫瘤免疫反應(yīng)，將“冷”免疫型腫瘤和無免疫反應(yīng)腫瘤轉(zhuǎn)化為“熱”免疫型腫瘤。此外，化療藥物一直在晚期惡性腫瘤的治療中占據(jù)著重要的地位?；熕幬锒嘁蕾噭┝繗麗盒阅[瘤細(xì)胞，劑量越大療效越好，但由于患者自身因素、藥物副作用等的限制，化療藥物的使用劑量往往有上限[14]。如果將誘導(dǎo)細(xì)胞焦亡與化療藥物協(xié)同使用，可以增加腫瘤細(xì)胞單次殺傷能力，減少化療藥物毒副作用，獲得更好的療效[14]。

總之，LUAD癌組織和癌旁組織之間存在著41個差異表達(dá)的焦亡基因，可用于區(qū)分LUAD癌組織和癌旁組織，表明細(xì)胞焦亡和LUAD存在緊密的聯(lián)系?；诮雇鱿嚓P(guān)預(yù)后模型獲得的風(fēng)險評分是預(yù)測LUAD預(yù)后的獨(dú)立危險因素，與腫瘤免疫有關(guān)。但是，本研究仍存在一定的局限，未來仍需要進(jìn)一步的實(shí)驗(yàn)來驗(yàn)證預(yù)后模型在LUAD病例中的臨床價值。