陳睿鵬，臧洪婧，孫 意

（中南大學(xué)湘雅二醫(yī)院病理科，湖南 長(zhǎng)沙 410011）

乳腺癌（breast cancer，BC）是女性最常見(jiàn)的惡性腫瘤，是一組高度異質(zhì)性的疾病[1]。隨著生活方式的改變，自20 世紀(jì)70 年代以來(lái)，乳腺癌在世界范圍內(nèi)的發(fā)病率顯著上升[2]。焦亡（pyroptosis）也稱細(xì)胞炎性死亡，定義為炎性介質(zhì)介導(dǎo)的程序性死亡[3]。細(xì)胞焦亡的特點(diǎn)是gasdermin 被啟動(dòng)形成細(xì)胞膜孔，這與細(xì)胞凋亡機(jī)制不同[4]。在促炎因子作用下，完整的細(xì)胞膜被破壞，最終導(dǎo)致細(xì)胞死亡[5]。研究表明[6]，焦亡在人類惡性腫瘤的發(fā)生發(fā)展中起不可或缺的作用。焦亡與癌癥的作用機(jī)制是復(fù)雜的[7，8]：一方面，介導(dǎo)腫瘤細(xì)胞焦亡以抑制其生長(zhǎng)；另一方面，焦亡為腫瘤的生長(zhǎng)提供了營(yíng)養(yǎng)和免疫環(huán)境支持[9，10]。已有研究發(fā)現(xiàn)[11-13]，焦亡相關(guān)基因與卵巢癌、肺腺癌、胃癌有關(guān)。這些研究表明，焦亡基因可作為一種新的基因標(biāo)記預(yù)測(cè)患者的生存預(yù)后。然而，焦亡基因在乳腺癌中的表達(dá)水平、與臨床預(yù)后特征和免疫浸潤(rùn)狀態(tài)的關(guān)系尚未明確。為此，本研究利用癌癥基因組圖譜（The Cancer Genome Atlas，TCGA）的基因表達(dá)資料構(gòu)建預(yù)后風(fēng)險(xiǎn)模型，擬評(píng)估乳腺癌患者的生存狀態(tài)，并使用基因表達(dá)綜合（Gene Expression Omnibus，GE0）數(shù)據(jù)庫(kù)進(jìn)行驗(yàn)證。

1 資料與方法

1.1 資料來(lái)源 TCGA-乳腺癌隊(duì)列的基因表達(dá)RNA測(cè)序數(shù)據(jù)及臨床信息從genome data Commons data Portal（https://portal.gdc.cancer.gov/）獲得。驗(yàn)證集GSE20685 和GSE42568 的數(shù)據(jù)來(lái)自GEO 數(shù)據(jù)門戶網(wǎng)站（https://www.ncbi.nlm.nih.gov/geo/）。在去除沒(méi)有臨床信息的數(shù)據(jù)后，本研究的訓(xùn)練集——TCGA 隊(duì)列納入了1045 個(gè)乳腺癌樣本和113 個(gè)相鄰正常樣本。驗(yàn)證集GEO 隊(duì)列納入了431 個(gè)乳腺癌樣本。此外，179 份正常人類乳腺組織樣本的基因表達(dá)RNA-seq 數(shù)據(jù)從GTEx 數(shù)據(jù)集（https://xenabrowser.net/datapages/）中獲得。TCGA 數(shù)據(jù)庫(kù)中的體細(xì)胞突變圖譜來(lái)源于genome Data Commons Data Portal（https://portal.gdc.cancer.gov/）。

1.2 焦亡相關(guān)差異表達(dá)基因及通路的關(guān)鍵基因 使用TCGA 數(shù)據(jù)進(jìn)行差異基因表達(dá)分析，發(fā)現(xiàn)乳腺癌樣本（1045 個(gè)）與癌旁正常組織樣本（113 個(gè)）的數(shù)量級(jí)存在差異，這可能導(dǎo)致統(tǒng)計(jì)誤差。因此，GTEX 數(shù)據(jù)集中正常乳腺組織樣本（179 個(gè)）也被納入研究，以確保結(jié)果的正確性。該網(wǎng)站（https://xenabrowser.net/datapages/）已對(duì)數(shù)據(jù)預(yù)處理，去除了批次效應(yīng)，可直接運(yùn)用。此外，在線數(shù)據(jù)庫(kù)STRING（v11.0,http://www.string-db.org/）用于可視化焦亡相關(guān)基因的蛋白質(zhì)相互作用網(wǎng)絡(luò)（PPI），PPI 分析所需的最小交互評(píng)分設(shè)置為0.4（默認(rèn)推薦值）[14]。將結(jié)果文件導(dǎo)入Cytoscape（https://cytoscape.org/）通過(guò)MNC 算法鑒定焦亡通路的關(guān)鍵基因（hubgene）[15，16]。

1.3 構(gòu)建與驗(yàn)證風(fēng)險(xiǎn)評(píng)分模型 使用Kaplan-Meier曲線和Log-rank 檢驗(yàn)從上述步驟得到的焦亡相關(guān)差異表達(dá)基因集中篩選出具有預(yù)后價(jià)值的基因（P＜0.05）。然后，基于具有顯著預(yù)后價(jià)值的焦亡相關(guān)基因進(jìn)一步分析。應(yīng)用最小絕對(duì)收縮和選擇操作數(shù)（LASSO）懲罰方法構(gòu)建預(yù)后模型。通過(guò)10 倍交叉驗(yàn)證確定最佳參數(shù)λ。根據(jù)“survminer”R 包計(jì)算的最佳截?cái)嘀?，將TCGA 訓(xùn)練集分為高、低風(fēng)險(xiǎn)組，通過(guò)“timeROC”軟件包繪制依賴于時(shí)間的受試者工作特征（ROC）曲線，評(píng)估模型在TCGA 和GEO 隊(duì)列中的特異性和敏感性。

1.4 風(fēng)險(xiǎn)評(píng)分與臨床分子特征及免疫浸潤(rùn)的相關(guān)性使用R 語(yǔ)言，采用Wilcoxon 秩和檢驗(yàn)（分類變量為2 個(gè)）及Kruskal-Wallis 檢驗(yàn)（分類變量大于2 個(gè)）探究焦亡風(fēng)險(xiǎn)評(píng)分與臨床分子特征的相關(guān)性。隨后，使用TIMER 數(shù)據(jù)庫(kù)（https://cistrome.shiny apps.io/timer/）闡明TCGA 隊(duì)列中高、低風(fēng)險(xiǎn)組的免疫浸潤(rùn)水平差異[17]。

1.5 基于風(fēng)險(xiǎn)亞組的差異基因進(jìn)行通路富集及基因組分析 在R 語(yǔ)言中通過(guò)Wilcoxon 秩和檢驗(yàn)評(píng)估高、低風(fēng)險(xiǎn)亞組的差異表達(dá)基因。閾值設(shè)為錯(cuò)誤發(fā)現(xiàn)率（FDR）＜0.05，樣本質(zhì)檢表達(dá)量的差異倍數(shù)（FC）≥2。隨后，利用WebGestalt（https://www.webgestalt.org/option.php）對(duì)焦亡高、低風(fēng)險(xiǎn)組的差異表達(dá)基因進(jìn)行功能和通路富集分析，包括基因本體論（GO）和京都基因與基因組百科全書（KEGG）[18]。將目標(biāo)基因輸入WebGestalt 網(wǎng)頁(yè)，應(yīng)用“ggplot2”R 包繪制結(jié)果。此外，應(yīng)用R 包通過(guò)“maftools”R 包處理以突變注釋格式（MAF）形式排序的體細(xì)胞突變譜，探究高、低風(fēng)險(xiǎn)亞組的基因組改變情況。

1.6 列線圖的構(gòu)建與評(píng)估 單因素及多因素Cox 回歸分析評(píng)估風(fēng)險(xiǎn)評(píng)分及臨床特征與生存預(yù)后的關(guān)系。整合風(fēng)險(xiǎn)評(píng)分及有預(yù)后價(jià)值的臨床特征構(gòu)建列線圖（nomogram）。應(yīng)用1、3 和5 年生存率的校準(zhǔn)曲線來(lái)評(píng)估列線圖與理想模型的偏差。構(gòu)建ROC 曲線，計(jì)算曲線下面積（AUC），以確定列線圖生存預(yù)后的預(yù)測(cè)能力。

1.7 統(tǒng)計(jì)學(xué)方法 所有數(shù)據(jù)通過(guò)R 語(yǔ)言（版本4.0.5）和R Bioconductor 包進(jìn)行處理。本研究的主要終點(diǎn)是總生存期（OS），即從最初診斷到最終死亡的時(shí)間間隔。用Kaplan-Meier 曲線和Log-rank 檢驗(yàn)比較風(fēng)險(xiǎn)評(píng)分組的總生存期，使用“rms”包構(gòu)造列線圖，C 指數(shù)用“survcomp”軟件包計(jì)算。除非特別說(shuō)明，所有統(tǒng)計(jì)檢驗(yàn)均為雙向檢驗(yàn)，P＜0.05 為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 焦亡基因的差異表達(dá)分析及關(guān)鍵基因 共收集35 個(gè)焦亡基因，檢測(cè)出32 個(gè)為差異表達(dá)基因（P＜0.05），見(jiàn)圖1A；其中18 個(gè)基因（CASP1、CASP4、CASP8、CASP9、ELANE、GPX4、GSDMB、GSDMD、GSDME、IL6、NLRP1、NLRP3、NOD1、PJVK、PLCG1、PRKACA、SCAF11 和TIRAP）在腫瘤組織中表達(dá)下調(diào)，14 個(gè)基因（AIM2、CASP3、CASP6、GSDMA、GSDMC、GZMA、IL18、IL1B、NLRC4、NLRP2、NLRP7、NOD2、PYCARD 和TNF）在腫瘤組織中表達(dá)上調(diào)。此外，35個(gè)基因集中PYCARD、AIM2、IL18、CASP1、CASP5、CASP8、NLRC4、IL6、NLRP3 和TNF 為焦亡通路的關(guān)鍵基因，見(jiàn)圖1B。

圖1 焦亡相關(guān)基因在乳腺癌中的表達(dá)水平及蛋白相互作用

2.2 焦亡風(fēng)險(xiǎn)評(píng)分模型構(gòu)建 從32 個(gè)焦亡相關(guān)差異表達(dá)基因中鑒定出5 個(gè)具有預(yù)后價(jià)值的基因：CASP9、GSDMC、GZMA、IL18 和PYCARD，均與乳腺癌患者的預(yù)后有關(guān)，見(jiàn)圖2。基于這些具有預(yù)后價(jià)值的基因，使用LASSO-Cox 回歸分析構(gòu)建評(píng)分模型（圖3A），通過(guò)10 倍交叉驗(yàn)證（圖3B），得到最小lambda 值=0.000 753 102 3，將相應(yīng)系數(shù)代入，得到計(jì)算公式：風(fēng)險(xiǎn)評(píng)分=0.107 146 60×GSDMC 表達(dá)值-0.207 295 50×CASP9 表達(dá)值-0.205 471 10×IL18表達(dá)值-0.038 234 58×PYCARD 表達(dá)值-0.092 932 92×GZMA 表達(dá)值。由此可知，CASP9、IL18、PYCARD、GZMA 表示上調(diào)與預(yù)后良好相關(guān)，而GSDMC 表達(dá)上調(diào)與預(yù)后不良相關(guān)。使用最佳截?cái)嘀担╟ut-off=0.056 399 17）將TCGA 隊(duì)列分為高、低風(fēng)險(xiǎn)組，結(jié)果顯示高、低風(fēng)險(xiǎn)組的KM 生存曲線比較，差異有統(tǒng)計(jì)學(xué)意義（P＜0.05），見(jiàn)圖3C；高低風(fēng)險(xiǎn)組的生存死亡狀況、人群分布趨勢(shì)以及構(gòu)成風(fēng)險(xiǎn)評(píng)分的基因表達(dá)熱圖見(jiàn)圖3D；1、3、5 年的ROC 曲線的AUC 分別為0.65、0.61、0.62，見(jiàn)圖3E，表明風(fēng)險(xiǎn)評(píng)分預(yù)測(cè)乳腺癌的預(yù)后有較高敏感性和準(zhǔn)確性。

圖2 有預(yù)后價(jià)值的焦亡基因KM 曲線

圖3 乳腺癌焦亡基因預(yù)后模型的構(gòu)建

2.3 外部數(shù)據(jù)集驗(yàn)證模型 以GSE20685（n=327）和GSE42568（n=104）作為外部驗(yàn)證數(shù)據(jù)集，使用TCGA 訓(xùn)練集的最佳臨界值（cut-off=0.056 399 17）將GEO 驗(yàn)證隊(duì)列分為高、低風(fēng)險(xiǎn)組，KM 生存曲線顯示，高風(fēng)險(xiǎn)組預(yù)后差于低風(fēng)險(xiǎn)組，見(jiàn)圖4A、圖4B，其中兩個(gè)GEO 隊(duì)列的生存和死亡狀態(tài)、人群分布和預(yù)后相關(guān)焦亡基因的表達(dá)水平見(jiàn)圖4C、圖4D。在驗(yàn)證集GSE20685 中，1、3、5 年的ROC 曲線的AUC 分別為0.77、0.68、0.69，見(jiàn)圖4E；在驗(yàn)證集GSE42568中，1、3、5 年的ROC 曲線的AUC 分別為0.62、0.66、0.72，見(jiàn)圖4F。

圖4 GSE20685 和GSE42568 隊(duì)列驗(yàn)證風(fēng)險(xiǎn)評(píng)分

2.4 風(fēng)險(xiǎn)評(píng)分與臨床特征、分子亞型及免疫浸潤(rùn)的相關(guān)性 TCGA 隊(duì)列中高、低風(fēng)險(xiǎn)組的患者分布趨勢(shì)見(jiàn)圖5A～圖5D；男性、導(dǎo)管癌、腫瘤大于2 cm、三陰分子亞型、高M(jìn)SI 和TMB 狀態(tài)的患者風(fēng)險(xiǎn)評(píng)分較高（P＜0.05），見(jiàn)圖5E～圖5J。免疫浸潤(rùn)結(jié)果顯示，低風(fēng)險(xiǎn)組活化CD4 T 細(xì)胞、活化CD8 T 細(xì)胞等效應(yīng)細(xì)胞浸潤(rùn)相對(duì)較高，見(jiàn)圖5K。

2.5 高低風(fēng)險(xiǎn)組富集途徑和基因組改變分析 在TCGA 隊(duì)列的高風(fēng)險(xiǎn)組中，562 個(gè)基因上調(diào)，210 個(gè)基因下調(diào)，見(jiàn)圖6。GO 分析顯示，高風(fēng)險(xiǎn)組上調(diào)基因富集于上皮細(xì)胞的生長(zhǎng)及角化過(guò)程（圖7A），高風(fēng)險(xiǎn)組的下調(diào)基因富集于個(gè)別信號(hào)通路及肌球蛋白的合成（圖7B）。KEGG 分析顯示，上調(diào)基因主要富集在炎癥介質(zhì)介導(dǎo)的信號(hào)通路（圖7C），下調(diào)基因主要富集于生物代謝反應(yīng)（圖7D）。高低風(fēng)險(xiǎn)組最頻繁突變的基因并不完全一致（圖7E、圖7F），這揭示了焦亡風(fēng)險(xiǎn)評(píng)分組基因組改變的差異。采用Fisher 檢驗(yàn)對(duì)高低風(fēng)險(xiǎn)組不同突變基因進(jìn)行探究，結(jié)果發(fā)現(xiàn)TP53處于首位（圖7G），表明TP53 與乳腺癌發(fā)生細(xì)胞焦亡有較高的相關(guān)性且提示預(yù)后不良。在高低風(fēng)險(xiǎn)組共發(fā)生和互排斥的突變中，TP53-PIK3CA 存在互斥突變現(xiàn)象（圖7H、圖7I），提示兩者的突變?cè)谌橄侔┲锌赡墚a(chǎn)生拮抗作用。

圖6 差異基因的火山圖

圖7 高低風(fēng)險(xiǎn)亞組的富集途徑和基因組改變分析

圖7 高低風(fēng)險(xiǎn)亞組的富集途徑和基因組改變分析（續(xù)）

2.6 構(gòu)建列線圖 單因素和多因素Cox 分析表明，焦亡風(fēng)險(xiǎn)評(píng)分可作為預(yù)測(cè)乳腺癌生存狀態(tài)的獨(dú)立預(yù)后指標(biāo)，見(jiàn)表1。整合風(fēng)險(xiǎn)評(píng)分和預(yù)后相關(guān)的臨床特征（年齡、遠(yuǎn)程轉(zhuǎn)移分期、TNM 分期）以構(gòu)建列線圖（圖8A）。時(shí)間依賴性ROC 曲線顯示，列線圖預(yù)測(cè)效果優(yōu)于單獨(dú)的風(fēng)險(xiǎn)評(píng)分（圖8B）；列線圖的KM 生存曲線結(jié)果顯著（圖8C）；校準(zhǔn)曲線表明，與理想模型相比，列線圖性能良好（C-index=0.7618，圖8D）。

表1 TCGA 訓(xùn)練集的單變量和多變量分析

表1（續(xù)）

圖8 基于TCGA 隊(duì)列構(gòu)建列線圖

圖8 基于TCGA 隊(duì)列構(gòu)建列線圖（續(xù)）

3 討論

細(xì)胞焦亡是一種獨(dú)特的程序性細(xì)胞死亡，不同人體組織和遺傳景觀下焦亡對(duì)癌癥有不同影響[7]。本研究從已發(fā)表的文獻(xiàn)中收集了35 個(gè)焦亡基因，其中大多數(shù)為差異表達(dá)基因。與正常組織相比，腫瘤組織中18 個(gè)焦亡基因表達(dá)減少，14 個(gè)基因表達(dá)增加。PYCARD、AIM2、IL18、CASP1、CASP5、CASP8、NLRC4、IL6、NLRP3 和TNF 是焦亡通路的關(guān)鍵基因。本研究采用LASSO Cox 回歸分析建立了一個(gè)包括5 個(gè)焦亡基因的風(fēng)險(xiǎn)評(píng)分，該評(píng)分在預(yù)測(cè)乳腺癌患者1、3、5 年總生存率方面具有較高的準(zhǔn)確性和實(shí)用性。結(jié)果表明，CASP9、IL18、PYCARD、GZMA 是預(yù)后的保護(hù)因素，GSDMC 是預(yù)后的危險(xiǎn)因素。

CASP9 是一種半胱氨酸蛋白酶，可觸發(fā)細(xì)胞凋亡[19，20]。有研究表明[21]，在人類乳腺癌中，通過(guò)上調(diào)CASP9 表達(dá)抑制miR-182-5p 來(lái)誘導(dǎo)腫瘤細(xì)胞凋亡并抑制其增殖。IL18 參與了多種惡性腫瘤的發(fā)生發(fā)展，文獻(xiàn)報(bào)道[22]其在胃腸道癌、乳腺癌和惡性黑色素瘤中表達(dá)上調(diào)，這與本研究的結(jié)果一致。IL18 具有促腫瘤和抗腫瘤的雙重作用，在結(jié)腸炎相關(guān)結(jié)直腸癌中可通過(guò)修復(fù)上皮屏障對(duì)機(jī)體起保護(hù)作用[23]；而在肝細(xì)胞癌中可驅(qū)動(dòng)癌細(xì)胞發(fā)生轉(zhuǎn)移致預(yù)后不良[24]。促凋亡基因PYCARD 在焦亡通路中起重要作用，目前尚缺少其對(duì)乳腺癌患者預(yù)后影響的報(bào)道，但有研究表明相較于正常組織，乳腺癌患者中PYCARD 高表達(dá)。這很可能與PYCARD 基因的甲基化水平有關(guān)[25，26]。毒素顆粒酶A（GZMA）是一種通過(guò)Caspase 途徑介導(dǎo)細(xì)胞凋亡的胰蛋白酶[27]。GZMA 在大多數(shù)胰腺癌（70%）中表達(dá)，在＜35%的乳腺癌、宮頸癌、肝癌、卵巢癌、前列腺癌、腎癌、胃癌、睪丸癌和尿路上皮癌中表達(dá)，在＜10%的淋巴瘤和黑色素瘤中表達(dá)。CTL/NK 細(xì)胞可通過(guò)表達(dá)上調(diào)GZMA 殺滅癌細(xì)胞[28]。GSDMC 屬于gasdermin 家族中一員，有報(bào)道稱GSDMC 參與乳腺癌細(xì)胞焦亡并在該病中高表達(dá)，同時(shí)與患者預(yù)后不良相關(guān)[29]。

腫瘤微環(huán)境在抗腫瘤分子治療中起至關(guān)重要的作用。同在卵巢癌及胃癌觀察到的現(xiàn)象類似[12，13]，本研究中低風(fēng)險(xiǎn)組的免疫細(xì)胞浸潤(rùn)程度高于高風(fēng)險(xiǎn)組，這進(jìn)一步證實(shí)了焦亡在腫瘤免疫微環(huán)境中發(fā)揮著重要作用，而高風(fēng)險(xiǎn)組患者較差的預(yù)后可能是由于抗腫瘤免疫水平降低所致。比較高低風(fēng)險(xiǎn)組的基因組變化，提示高風(fēng)險(xiǎn)組與更具有侵略性的分子改變（如TP53 突變）顯著相關(guān)，同時(shí)本研究觀察到一種特殊的TP53-PIK3CA 互斥突變現(xiàn)象，這些基因組的改變與乳腺癌患者預(yù)后不良高度相關(guān)。本研究通過(guò)單變量和多變量Cox 分析發(fā)現(xiàn)，結(jié)合焦亡風(fēng)險(xiǎn)評(píng)分、年齡、遠(yuǎn)程轉(zhuǎn)移分期、TNM 分期的列線圖具有優(yōu)秀的預(yù)測(cè)能力，可作為一項(xiàng)有意義的個(gè)性化醫(yī)療指標(biāo)。

綜上所述，本研究構(gòu)建了焦亡相關(guān)基因的乳腺癌預(yù)后模型，該模型有利于臨床醫(yī)生制訂個(gè)性化的治療方案，提高乳腺癌患者的遠(yuǎn)期生存率。