肺腺癌焦亡相關特征模型的構建和驗證

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

（浙江中醫藥大學附屬溫州中醫院，浙江 溫州 325000）

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

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

1 材料和方法

1.1 數據收集

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

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

1.2 差異表達基因分析

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

1.3 差異表達基因的共識聚類分析

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

1.4 構建基于LUAD細胞焦亡基因的預后風險模型

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

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

1.5 風險評分和LUAD患者生存之間的關系

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

1.6 生物學功能分析

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

1.7 免疫浸潤分析

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

1.8 統計學分析

采用R 4.0.3軟件進行統計分析。

2 結果

2.1 差異表達的焦亡基因

共篩選出41個差異表達的焦亡基因，其中表達下調10個（IL6、NLRC4、CASP5、IL1B、NLRP3、IL1A、CASP1、IRF1、CHMP3、NLRP1）、表達上調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分析結果表明（相互作用分數為0.9），CHMP3、CHMP4C、CHMP2A、CHMP7、CHMP4B、CHMP2B、CHMP6、CHMP4A、CASP6、CASP3、CASP8是核心基因，見圖1（b）。以閾值r=0.5為臨界點繪制相關網絡圖，見圖1（c）。

2.2 基于差異表達基因的腫瘤分型

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

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

2.3 基于LUAD焦亡基因的預后模型構建及其臨床價值

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

圖3 預后模型對TCGA隊列的預后評估效能

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

圖4 預后模型對GEO隊列的預后評估效能

2.4 預后風險模型與LUAD患者預后的關系

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

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

2.5 TCGA隊列高風險組和低風險組差異表達基因的生物學功能分析

KEGG富集分析和GO富集分析結果顯示，TCGA隊列高風險組和低風險組差異表達基因主要與趨化因子介導腫瘤相關信號通路、免疫反應、細胞膜功能相關。見圖6。

圖6 TCGA隊列高風險組、低風險組差異表達基因的功能分析和ssGSEA結果

采用ssGSEA量化TCGA隊列LUAD患者的免疫相關通路。與高風險組比較，低風險組B淋巴細胞、抗原遞呈細胞、T淋巴細胞、人類白細胞抗原、中性粒細胞、腫瘤浸潤淋巴細胞、調節性T細胞的比例更高。見圖6。

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

3 討論

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

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

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

近年來，相關研究構建了很多LUAD的風險預后模型，為LUAD的診斷和預后提供了新的思路。本研究證實了焦亡基因對LUAD的診斷和預后具有重要價值，通過LASSO-Cox回歸分析建立風險預后模型，采用傳統方法證明風險預后模型的有效性。細胞焦亡會激活機體免疫反應，發揮抗癌的作用，設計誘導細胞焦亡的治療策略可以加強抗腫瘤免疫反應，將“冷”免疫型腫瘤和無免疫反應腫瘤轉化為“熱”免疫型腫瘤。此外，化療藥物一直在晚期惡性腫瘤的治療中占據著重要的地位。化療藥物多依賴劑量殺傷惡性腫瘤細胞，劑量越大療效越好，但由于患者自身因素、藥物副作用等的限制，化療藥物的使用劑量往往有上限[14]。如果將誘導細胞焦亡與化療藥物協同使用，可以增加腫瘤細胞單次殺傷能力，減少化療藥物毒副作用，獲得更好的療效[14]。

總之，LUAD癌組織和癌旁組織之間存在著41個差異表達的焦亡基因，可用于區分LUAD癌組織和癌旁組織，表明細胞焦亡和LUAD存在緊密的聯系。基于焦亡相關預后模型獲得的風險評分是預測LUAD預后的獨立危險因素，與腫瘤免疫有關。但是，本研究仍存在一定的局限，未來仍需要進一步的實驗來驗證預后模型在LUAD病例中的臨床價值。