PDF

J Cancer 2024; 15(2):428-443. doi:10.7150/jca.90113 This issue Cite

Research Paper

Identification of a PANoptosis-related Gene Signature for Predicting the Prognosis, Tumor Microenvironment and Therapy Response in Breast Cancer

Tianshui Yu, Weilun Cheng, Jiarui Zhang, Ting Wang, Yansong Liu, Yunqiang Duan, Anbang Hu, Jianyuan Feng, Mingcui Li, Yanling Li, Hanyu Zhang, Zhiyuan Rong, Yuhang Shang, Suborna S. Shakila, Fei Ma , Baoliang Guo

Department of General Surgery, The Second Affiliated Hospital of Harbin Medical University, Harbin 150081, China.

✉ Corresponding authors: Baoliang Guo (baoliangguo2020edu.cn); Fei Ma (wafsfdcom).

Abstract

Breast cancer (BC) is the most prevalent malignancy among women worldwide. Mounting evidence suggests that PANoptosis participates in cancer development and therapy. However, the role of PANoptosis in BC remains unclear. In this study, we identified ten PANoptosis-related genes using Cox regression analysis, random forest (RF) algorithm and least absolute shrinkage and selection operator (LASSO) algorithm. A PANoptosis-related score (PRS) was calculated based on the coefficient of LASSO. Notably, we divided the patients into high- and low-risk groups according to the PRS and revealed a negative correlation between PRS and overall survival. Next, a nomogram model was constructed and validated to improve the clinical application of PRS. Functional enrichment analyses and the Bayesian network demonstrated that differentially expressed genes between high- and low-risk groups were mainly enriched in immune-related pathways. Besides, we found significant differences in tumor mutation burden and tumor immune microenvironment between patients in these two groups using bulk-RNA and single-cell RNA sequencing data. Furthermore, charged multivesicular body protein 2B (CHMP2B) was identified as the hub gene by combining LASSO, weighted gene co-expression network analysis, RF and eXtreme Gradient Boosting. Importantly, using immunohistochemistry analysis based on our tissue microarray, we found that CHMP2B was highly expressed in tumor tissue, and CD4 and CD8 were more likely to be positive in the CHMP2B-negative group. Survival analyses revealed that CHMP2B adversely impacted the survival of BC patients. In conclusion, we not only constructed a highly accurate predictive model based on PRS, but also revealed the importance of PANoptosis-related gene signature in the modulation of the tumor microenvironment and drug sensitivity in BC.

Keywords: breast cancer, PANoptosis, prognostic model, tumor microenvironment, drug sensitivity, machine learning

Introduction

Breast cancer (BC) is the most prevalent primary malignancy among women worldwide, accounting for 1 in 4 cancer cases and 1 in 6 cancer deaths of all cancers [1]. Currently, BC treatment approaches mainly include surgery, chemotherapy, targeted therapy, endocrine therapy and radiotherapy [2]. Considering the heterogeneity and complexity of BC, it is crucial to screen high-risk patients and implement more appropriate treatment strategies based on pathological characteristics [3]. With the development of the DNA microarray and next-generation sequencing (NGS) over the past decades, individual gene signatures could provide alternative information to predict treatment sensitivity and prognosis of BC patients in addition to clinicopathological features [4-6].

Resistance to cell death is a hallmark of cancer [7]. Cell death can be divided into accidental cell death (ACD) and regulated cell death (RCD), and the latter is generally referred to as programmed cell death (PCD) under physiological conditions [8]. Apoptosis, pyroptosis and necroptosis are three critical PCD pathways characterized by particular molecular and genetic features [8]. Many studies have already explicitly elucidated the importance of these pathways in the carcinogenesis and treatment of BC [9-11]. Recently, based on the extensive cross-talk between these PCD pathways, a conceptualization of an integrated cell death modality called “PANoptosis” was formed [12]. PANoptosis is an inflammatory PCD pathway activated by the simultaneous involvement of components from pyroptosis, apoptosis and necroptosis [13]. It could not be characterized by any of these pathways alone, although it has the essential features of each of them [14]. PANoptosis has been proven to participate in regulating the tumorigenesis of colorectal cancer and the immunotherapy response of gastric cancer [15, 16]. Besides, He et al. found that PANoptosis was also related to the survival of BC patients [17]. However, more comprehensive research is still needed to better unravel its role in modulating BC progression and therapeutic response, as well as the ability to predict the prognosis of BC patients in combination with other clinical information.

In this study, we constructed a nomogram based on PANoptosis-related genes and clinical features to predict the prognosis of BC patients for the first time, not only identifying the high-risk patients, but also helping to implement more appropriate treatment for certain patients.

Materials and Methods

Data source

The gene expression profile and related clinicopathological data of TCGA-BRCA were downloaded from UCSC Xena and utilized as the training cohort in this study [18]. The METABRIC breast cancer data was downloaded from cBioPortal for external validation [19, 20]. All of the patients enrolled had primary breast cancer and M0 tumor stage. The baseline characteristics are shown in Table S1 and compared using the chi-square test or Fisher's exact test. PANoptosis-related genes were created by combining the gene lists of pyroptosis, apoptosis and necroptosis, which were collected from the Molecular Signatures Database (MsigDB) and literature review [21]. Single-cell RNA sequencing (scRNA-seq) data was extracted from the GSE161529 dataset in the Gene Expression Omnibus (GEO) database [22]. One triple-negative breast cancer (TNBC), one human epidermal growth factor receptor 2-positive (HER2+) and one estrogen receptor-positive (ER+) cases were selected for the study. The workflow of this study is shown in Figure S1.

Construction of PANoptosis-related risk score and nomogram

We used the Cox proportional hazard regression model to assess the association between the expression of each PANoptosis-related gene and the overall survival (OS) of patients in the TCGA-BRCA cohort. Afterward, the random forest (RF) algorithm was conducted to further screen featured genes using the randomForest R package [23]. Genes with MeanDecreaseGini > 10 were selected for further analysis. The Least absolute shrinkage and selection operator (LASSO) algorithm was then applied to ensure the model's simplicity and minimize overfitting in the model training process [24]. The risk score was constructed by using the regression coefficients derived from LASSO Cox regression analysis:

PRS means the PANoptosis-related score. Coef_i and Exp_i represent the coefficient and expression level of the corresponding gene, respectively. Based on the median risk score, all samples were separated into two categories: high-risk group (PRS > median value) and low-risk group (PRS < median value). Typical immunohistochemical (IHC) images of all model genes were downloaded from the Human Protein Atlas (HPA) database. The Kaplan-Meier (KM) analysis was used to compare the differences in OS between these two groups. The time-dependent receiver operating characteristic (ROC) curve was generated to evaluate the prognostic performance of the PRS. Combined with the clinicopathological features that were statistically significant in univariate (p-value < 0.10) and multivariate (p-value < 0.05) Cox regression analysis, we subsequently used the “rms” package to predict 1-, 3-, 5-, and 10-year OS by constructing a nomogram. The nomogram's discrimination performance was quantitatively assessed by the area under curve (AUC) of the ROC curve, and the calibration performance was evaluated by the calibration curve.

Functional enrichment analysis

Differentially expressed genes (DEGs) between high- and low-risk groups were identified using the “DESeq2” package in R software [25]. Functional enrichment analyses based on Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene Ontology (GO) and Reactome databases were carried out between these two groups, and the results were evaluated by the R package “clusterProfiler” [26] and “ReactomePA” [27]. The Bayesian network was applied using the “CBNplot” R package to understand how GO and Reactome pathways interact with one another [28].

Somatic mutation analysis

Somatic mutation data retrieved from The Cancer Genome Atlas (TCGA) were analyzed using the R package “maftools” [29]. The tumor mutational burden (TMB) of each patient between high- and low-risk groups was calculated and compared using the Wilcoxon rank-sum test.

Immunogenomic landscape analysis

We used CIBERSORT to predict the proportions of 22 types of tumor-infiltrating immune cells (TIICs) in each sample. Single sample Gene Set Enrichment Analysis (ssGSEA) was applied to predict the abundance of 28 TIICs in individual tissue samples [30]. Furthermore, we used the R package “estimate” to calculate each sample's ESTIMATE score. The Immunophenoscore (IPS) scores, which are calculated based on representative cell-type gene expression z-scores, were collected from the Cancer Immunome Atlas (TCIA) [31]. The tumor immune exclusion score was calculated by the Tumor Immune Dysfunction and Exclusion (TIDE) [32, 33]. Moreover, we compared the expression of genes from the human leukocyte antigen (HLA) family and several critical co-stimulators between high- and low-risk groups using the Wilcoxon rank-sum test.

Identification of PANoptosis-associated single cells

Seurat (v4.3.0) was utilized for quality control and processing of scRNA-seq data [34]. Cell types were annotated based on the CellMarker 2.0 database and relative literature [35, 36]. Single-Cell Map of Diverse Immune Phenotypes in the Breast Tumor Microenvironment (Scissor) package was used to identify bulk phenotype-associated cell subpopulations [37]. In this study, the input data of the Scissor pipeline included TCGA-BRCA bulk data, GSE161529 scRNA-seq data and PANoptosis-related risk group data obtained from LASSO analysis. Based on the signs of the estimated regression coefficients, the single cells were finally grouped into Scissor positive (Scissor+) and Scissor negative (Scissor-) cells, corresponding to the high-risk group and low-risk group in TCGA-BRCA bulk data, respectively.

Therapeutic efficacy estimation

We utilized the R package “oncoPredict” to calculate the half-maximal inhibitory concentration (IC50) based on the GDSC2 dataset of the Genomics of Drug Sensitivity in Cancer (GDSC) database [38]. Besides, we assessed the expression of target genes associated with drugs exhibiting different sensitivities between high- and low-risk groups using the DrugBank database [39].

Screening for the hub gene

Weighted gene co-expression network analysis (WGCNA) was applied to identify hub genes among ten PANoptosis genes. We used the R package “WGCNA” to construct the co-expression network of genes from the TCGA-BRCA dataset [40]. According to their expression patterns, genes were classified into four modules (brown, blue, turquoise and grey). Subsequently, the correlation between modules and PRS was calculated. The brown module, which showed the highest correlation with PRS, was selected for further analysis. Besides, we applied eXtreme Gradient Boosting (XGBoost) to further calculate the importance of each gene in this ten-gene signature using the R package “XGBoost” [41]. Finally, the hub gene was characterized by combining the results from WGCNA, LASSO, RF and XGBoost.

Validation of the hub gene by IHC

The tissue microarray (TMA) was made by the production company (SHANGHAI OUTDO BIOTECH CO., LTD.) using wax blocks of 105 breast cancer and 41 normal breast tissue. The tissues were collected from patients with primary breast cancer and without metastases at diagnosis in the Second Affiliated Hospital of Harbin Medical University (Table S1). After de-paraffinization and rehydration, tissue sections were incubated in antigen retrieval buffer and heated in a steamer above 97 °C for 20 min. IHC staining was conducted utilizing Ventana Discovery XT Automated Slide Stainer (Ventana Medical Systems, Inc., Tucson, AZ, USA). The automated Discovery XT system was employed to carry out deparaffinization, antigen retrieval, blocking, DAB detection, counterstaining, post-counterstaining and slide cleaning. Charged multivesicular body protein 2B (CHMP2B) antibody (Abclonal, A19244), CD4 antibody (MXB biotechnologies, RMA-0620) and CD8 antibody (MXB biotechnologies, MAB-1031) were respectively applied to tissue sections overnight in a humidity chamber at a dilution of 1:200 at 4 °C. After washing in TBS, the antigen-antibody binding was detected using the Envision+ system and DAB+ chromogen (DAKO). Tissue sections were briefly immersed in hematoxylin for counterstaining, washed with water and covered with coverslips. Subsequently, slides were independently evaluated by two pathologists under a multi-headed microscope in a blinded manner. The staining intensity level was scored from 0 to 3 (no staining, light brown, brown, and tan). The staining extent was rated from 0 to 4 based on the percentage of positive cells (0-5%, 5-25%, 26-50%, 51-75%, and 76-100%). Finally, the IHC score was defined as the product of staining intensity and extent scores, ranging from 0 to 12. Scores 0-4 were grouped as negative and 5-12 as positive. Clinical characteristics were compared using Chi-square tests. Continuous variables were converted to categorical using X-tile software (Version 3.6.1) [42]. We then used the KM curve to estimate survival outcomes, and the Cox proportional hazards model was performed to fit clinical characteristics for OS. This research was approved by the Ethics Committee of Harbin Medical University (reference number KY2021-058).

Statistical analysis

All statistical analyses were implemented using R v4.2.2 and SPSS 26.0. A p-value < 0.05 was considered significant unless stated otherwise.

Results

Identification of a ten-gene PANoptosis-related signature in breast cancer

Sixty-five PANoptosis-related genes (Table S2) were input into univariate Cox regression analysis in the TCGA-BRCA training set, and thirteen genes significantly affected OS were obtained (Figure 1A). Subsequently, we applied the RF analysis and found that the MeanDecreaseGini for each of these genes was higher than 10 (Figure 1B). Then, the LASSO regression analysis was performed on these genes, and ten PANoptosis-related genes with non-zero coefficients were selected for further model construction (Figure 1C, D). Finally, based on the coefficients derived from the LASSO analysis, the formula for the risk score was as follows: PRS = (-0.174034262515715)*CHMP6 + (-0.0701780204898664)*TP63 + (-0.449396144853978)*IRF2 + (-0.121272009866694)*CASP7 + (-0.146272576423145)*CHMP4A + (-0.0749366454708963)*IL18 + 0.268668804297125*AIFM1 + 0.114110313114337*CHMP2B + (-0.0789991578059556)*GZMB + 0.179084139171232*CHMP3. To better unearth the potential mechanism of PANoptosis-related genes, we divided patients into high-risk and low-risk groups based on their PRS. All ten PANoptosis-related genes used for model construction were differentially expressed between the high- and low-risk groups (Figure 1E). Besides, we obtained typical IHC images of these genes from the HPA database (Figure 1F). Notably, the KM analysis further proved the prognostic value of PRS (Figure 1G). For external validation, the result of METABRIC data was consistent with the TCGA cohort (Figure 1H).

Construction and validation of a clinical nomogram

To develop a nomogram based on PRS and clinical features, we first applied time-dependent ROC to evaluate the prediction efficiency of PRS. As shown in Figure 2A, the AUC at 1-year, 3-year, 5-year and 10-year OS was 0.733, 0.694, 0.678 and 0.615, respectively. Next, we evaluated the prognostic effect of PRS, age, estrogen receptor (ER) status, T stage and lymph node status in univariate Cox regression analysis, and factors with p-value < 0.1 were enrolled in multivariate Cox regression analysis (Figure 2B). According to the results of multivariate Cox regression analysis, we selected age, ER status, lymph node status and PRS to develop the final model (Figure 2C). Given the complexity of the risk score formula, we constructed a nomogram to predict 1-, 3-, 5- and 10-year OS of BC patients (Figure 2D). In our model, higher OS rates were associated with a lower risk score, younger age, positive ER status and negative lymph node metastasis. The AUC of the nomogram was 0.857, 0.750, 0.721, and 0.804 for 1-, 3-, 5- and 10-year survival, respectively (Figure 2E). The calibration curve revealed that the predicted curve of the nomogram was nearly identical to the actual curve of BC patients, suggesting the close relationship between predicted survival rates and actual rates at 1, 3, 5 and 10 years (Figure 2F). In addition, the ROC and calibration curves of the external validation cohort further confirmed the significant predictive value of our model (Figure 2G, H). Taken together, these results demonstrated the outstanding performance of our nomogram compared with previously established prognostic models [43-45].

Identification of PANoptosis-related signaling pathways

We found 843 DEGs, including 275 up-regulated genes and 568 down-regulated genes with |log₂FC| > 1 and adjusted p-value < 0.05 between high- and low-risk groups (Figure 3A). In functional enrichment analysis, DEGs were mainly enriched in pathways correlated to immunity, including cytokine-cytokine receptor interaction, natural killer cell mediated cytotoxicity, T cell receptor signaling pathway, and Th1, Th2 and Th17 cell differentiation in KEGG analysis (Figure 3B, Table S3).

 Figure 1 

Identification of a ten-gene PANoptosis-related signature. (A) Forest plot of statistically significant genes in the univariate Cox regression analysis; (B) The MeanDecreaseGini index for the thirteen PANoptosis-related genes in the random forest (RF) analysis. (C) Coefficient profiles of PANoptosis-related genes in the least absolute shrinkage and selection operator (LASSO); (D) identification of the best parameter (lambda) in LASSO; (E) Violin plot of the expression of ten PANoptosis-related genes in high and low PRS groups; (F) typical immunohistochemical (IHC) images of ten PANoptosis-related genes in normal and breast tumor tissues from the HPA database; (G) Kaplan-Meier (KM) analysis based on PANoptosis-related score (PRS) in the training cohort; (H) KM analysis based on PRS in the validation cohort. HR: hazard ratio; CI: confidence interval; high: high-risk group; low: low-risk group; N: normal; T: tumor; ****: p < 0.0001.

GO analysis also revealed that DEGs were enriched in gene sets associated with T cell activation, lymphocyte mediated immunity, leukocyte cell-cell adhesion and immune receptor activity (Figure 3C, Table S4). Besides, enrichment based on the Reactome database demonstrated that DEGs were related to G protein-coupled receptor (GPCR) ligand binding, signaling by interleukins, T cell antigen receptor (TCR) signaling and PD-1 signaling (Figure 3D, Table S5). Furthermore, the BN showed that most pathways were intensely correlated (Figure 3E, F).

Tumor mutation analysis

Considering that genetic mutations are crucial factors in tumorigenesis, we estimated the situation of somatic mutation between two groups. The somatic mutation rate of the high-risk group was 86.89% (411 of 473 samples), primarily the missense mutation, and TP53 showed the highest frequency of mutations (40%) (Figure 4A). The mutation rate in the low-risk group was 87.86% (398 of 453 samples), mainly the missense mutation, and PIK3CA demonstrated the highest frequency of mutations (39%) (Figure 4B). In addition, TMB quantification analysis revealed a higher TMB in the high-risk group (Figure 4C-E).

Comprehensive analysis of tumor microenvironment

To delineate the immune status of the tumor microenvironment (TME), we performed the CIBERSORT algorithm to calculate the proportion of immune cells in each group (Figure 5A).

 Figure 2 

Construction and validation of a clinical nomogram. (A) Time-dependent receiver operating characteristic (ROC) curves of the PANoptosis-related score (PRS) in predicting the 1-, 3-, 5- and 10-year overall survival (OS); (B) Univariate Cox regression analysis of PRS and clinical characteristics; (C) Multivariate Cox regression analysis of PRS and clinical characteristics; (D) Nomogram based on the PRS, age, ER status and lymph node status; (E) Time-dependent ROC curves of the nomogram in predicting the 1-, 3-, 5- and 10-year OS in the training cohort; (F) Calibration curves of the nomogram in predicting the 1-, 3-, 5- and 10-year OS in the training cohort; (G) Time-dependent ROC curves of the nomogram in predicting the 1-, 3-, 5- and 10-year OS in the validation cohort; (H) Calibration curves of the nomogram in predicting the 1-, 3-, 5- and 10-year OS in the validation cohort. AUC, area under curve; HR: hazard ratio; CI: confidence interval; ER, estrogen receptor; OS, overall survival.

 Figure 3 

Identification of PANoptosis-related signaling pathways. (A) Volcano plot of DEGs based on PRS; (B) Dot plot of Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis; (C) Dot plot of Gene Ontology (GO) enrichment analysis; (D) Dot plot of Reactome enrichment analysis; (E) Bayesian network (BN) based on GO analysis; (F) Bayesian network (BN) based on Reactome analysis.

 Figure 4 

Tumor mutation analysis. (A) Waterfall plot of somatic mutation features established with high PANoptosis-related score (PRS); (B) Waterfall plot of somatic mutation features established with low PRS; (C) tumor mutation burden (TMB) in the high-risk group; (D) TMB in the low-risk group; (E) Violin plot of TMB in the high- and low-risk groups. high: high-risk group; low: low-risk group; ***, p < 0.001.

As expected, the high-risk group was associated with saliently fewer CD8+ T cells, M1 macrophages and more M2 macrophages than the low-risk group. ssGSEA also demonstrated a lower expression of most TIICs, including activated B cell, CD4+ T cell, CD8+ T cell and dendritic cell in the high-risk group (Figure 5B). Furthermore, the immune score calculated by the ESTIMATE algorithm was statistically lower in the high-risk group (Figure 5C). The results obtained from TCIA illustrated that the relative probabilities of responding to immunotherapy in the high-risk group were lower than those in the low-risk group, regardless of CTLA4 and PD1 status (Figure 5D). The high-risk group also showed a higher tumor immune exclusion score, indicating an immune suppressive microenvironment (Figure 5E). Moreover, most of the HLA genes (Figure 5F) and immune checkpoints (Figure 5G) were considerably higher expressed in the low-risk group.

PANoptosis-associated single cell analysis

We utilized the standard Seurat pipeline to explore the scRNA-seq data and got seven cell types after dimension reduction, cell clustering and annotation (Figure S2A, B, Figure 6A). Subsequently, the Scissor algorithm was applied to identify the PANoptosis-associated cell subpopulations. In this analysis, 884 Scissor+ cells (corresponding to the high-risk group) and 1021 Scissor- cells (corresponding to the low-risk group) were classified by the PRS of bulk samples (Figure 6B). Notably, we found that proportional fractions of T cell, B cell, mast cell and endothelial cell in the Scissor- group were higher than Scissor+ group, while epithelial cell, macrophage and fibroblast were lower (Figure 6C). Then, we extracted T cells, B cells and macrophages to further evaluate the relationship between PANoptosis and tumor immune microenvironment. In the subgroup analysis, immune cells were classified into ten types (Figure S2C, D, Figure 6D). Similarly, we got 458 Scissor+ cells and 520 Scissor- cells using the Scissor algorithm (Figure 6E). Concordantly, the proportion of cell types demonstrated that proportional fractions of effector memory CD8+ T cell, central memory CD4+ T cell, M1 macrophage, monocyte, B cell and plasma cell in the Scissor- group were higher than Scissor+ group, while T follicular helper cell, M2 macrophage, regulatory T cell and plasmacytoid dendritic cell were lower (Figure 6F). These results were consistent with our immune infiltration analysis using bulk data, further validating the positive correlation between higher PRS (high-risk group) and tumor-promoting immune microenvironment.

 Figure 5 

Comprehensive analysis of tumor immune microenvironment. (A) Boxplot of immune cell proportion in the high- and low-risk groups calculated by CIBERSORT algorithm; (B) Boxplot of immune cell expression in the high- and low-risk groups calculated by ssGSEA algorithm; (C) Violin plot of ESTIMATE score in the high- and low-risk groups; (D) Violin plot of IPS in the high- and low-risk groups; (E) Violin plot of TIDE exclusion score in the high- and low-risk groups; (F) Violin plot of the expression levels of HLA molecules; (G) Violin plot of the expression levels of co-stimulators. high: high-risk group; low: low-risk group; IPS, Immunophenoscore; ns, non-significant; *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

 Figure 6 

PANoptosis-associated TME analysis at single-cell resolution. (A) UMAP visualization of cell-type-specific annotation. (B) UMAP visualization of Scissor+ and Scissor- cells. (C) Proportional fractions of identified cell types across Scissor+/- condition. (D) UMAP visualization of cell-type-specific annotation among extracted immune cells. (E) UMAP visualization of Scissor+ and Scissor- cells among extracted immune cells. (F) Proportional fractions of identified cell types across Scissor+/- condition among extracted immune cells.

Drug sensitivity analysis

We explored the association between PRS and anticancer drug sensitivity in the GDSC database. Strikingly, we found some commonly used anticancer chemotherapeutic drug agents, including docetaxel (Figure 7A), epirubicin (Figure 7B), paclitaxel (Figure 7C), vinblastine (Figure 7D), vincristine (Figure 7E) and vinorelbine (Figure 7F), had a lower IC50 in the low-risk group, indicating their better responses in treating the low-PRS patients. We further assessed the expressions of target genes associated with the drugs mentioned above. Target genes from the DrugBank database, included JUN, MAP2, NR1I2, TOP2A, TUBB, TUBE1 and TUBG1, were differentially expressed between the high-risk and low-risk groups (Figure 7G). These results suggested that PRS could identify more suitable patients for appropriate anticancer drug therapy.

Hub gene screening and experimental verification

To identify the most crucial hub gene among ten model genes, we first initiated the WGCNA algorithm. A weighted gene co-expression network was established utilizing the TCGA-BRCA dataset (Figure S3A). The scale-free network was constructed by setting the soft threshold to 6 (Figure S3B). Moreover, we created an adjacency matrix and converted it into a Topological Overlap Matrix (TOM). Four gene modules were identified according to the TOM, namely blue (3465), brown (2684), grey (2599), and turquoise (10775) modules. We found that the brown module exhibited the strongest correlation with PRS (Figure S3C). Next, we extracted the overlapping genes between the brown module and the ten model genes, resulting in the discovery of four genes - CHMP2B, CHMP3, CHMP4A, and CHMP6. Finally, after incorporating the findings from the RF (Figure 1B) and XGBoost (Figure S3D) analysis, we selected CHMP2B as the hub gene for further experimental verification.

We then performed IHC analyses based on a TMA made by 105 BC and 41 normal breast tissue. In this cohort, 57.1% of BC patients were CHMP2B positive, while only 9.8% of healthy individuals highly expressed CHMP2B (Figure 8A, Table 1).

BC patients were then assigned to the CHMP2B-positive or CHMP2B-negative group according to their IHC scores. Unsurprisingly, we found that CD4 and CD8 were more likely to be positive in the CHMP2B-negative group, which was consistent with our former bioinformatical findings (Table 2).

 Figure 7 

Drug sensitivity analysis. (A-F) Sensitivity analysis for Docetaxel (A), Epirubicin (B), Paclitaxel (C), Vinblastine (D), Vincristine (E) and Vinorelbine (F) in patients at low and high risk; (G) Boxplot of target genes of differentially sensitive drugs. IC50, half-maximal inhibitory concentration; ns, non-significant; high: high-risk group; low: low-risk group; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

 Figure 8 

Experimental validation using IHC analysis. (A) Pie chart of CHMP2B expression in normal and breast tumor tissue. (B) Representative images of CHMP2B expression in normal (left panel) and tumor (right panel) tissue. (C) KM survival analysis between patients in CHMP2B-negative and CHMP2B-positive groups. IHC, immunohistochemistry; KM, Kaplan-Meier.

Typical photos of IHC in our cohort are shown in Figure 8B and Figure S4. Notably, the KM plot illustrated that patients in the CHMP2B-negative group shared better prognoses than those in the CHMP2B-positive group (Figure 8C). The univariate and multivariate Cox analysis further confirmed that CHMP2B was a detrimental prognostic factor in breast cancer (Table 3).

Discussion

Several types of PCD pathways, including apoptosis, necroptosis and pyroptosis, have been found to play an essential regulatory role in cancer development [46]. Apoptosis is generally divided into an intrinsic pathway and an extrinsic pathway, both mediated by the activation of initiator and executioner caspases [47]. Activation of necroptosis is mediated by the RIP kinases RIPK1 and RIPK3, together with the pore-forming pseudokinase MLKL, downstream of necrosome formation [48, 49]. Pyroptosis is activated through the inflammatory caspases (caspase-1 and caspase-4/5) and executed by the gasdermin family members (GSDMD and GSDME) [50]. PANoptosis is an integrated system in which any of the three PCD pathways can compensate for one another and work together at different times based on the context of the stimulus provided [14]. It is regulated by the PANoptosome complex, a molecular scaffold for engaging pivotal pyroptotic, apoptotic, and necroptotic machinery [12]. In oncology research, PANoptosis has been found to regulate the tumorigenesis of colorectal cancer and the immunotherapy response of gastric cancer [15, 16]. Additionally, PANoptosis-related genes could also predict the survival of BC patients and modulate tumor immune microenvironment [17]. Moreover, Maitra et al. reviewed the molecules in the PANoptosis pathway and proposed the hypothesis that PD-1/PD-L1-targeted inhibitor might play its role through PANoptosis pathway, but the specific mechanism of PANoptosis in breast cancer still needs further exploration [51].

Despite the similar research topic [17], this study screened PANoptosis-related genes using the combination of Cox regression, random forest and LASSO algorithm and built a highly accurate nomogram based on PANoptosis gene signature and clinicopathological information for the first time. Among genes in our model, a previous study has shown that tumor protein p63 (TP63) could act both as a tumor suppressor and an oncogene in breast cancer, depending on the cellular context [52]. Moreover, Interleukin 18 (IL-18), a pro-inflammatory cytokine, is up-regulated on tumor-infiltrating lymphocytes, suggesting that IL-18 therapy could enhance anti-tumor immunity [53].

Recently, several models have been developed to predict the prognosis of BC patients. Tian et al. have built a prognostic model using basement membrane-related genes [43]. Cui et al. have also established a nomogram based on nicotinamide metabolism-related signature [44], while Li's model was constructed utilizing m6A-related genes [45]. However, none of these models perform better than our PANoptosis-based nomogram (assessed using AUC), further proving the vital value of our research.

We not only constructed a PRS-based prognostic model, but also conducted a series of comprehensive analyses between patients in the high- and low-risk groups. We found that DEGs were enriched in several pathways associated with immunity in the functional enrichment analysis. Moreover, BN demonstrated an intense correlation among these pathways. These results suggested that different survival outcomes between patients in high- and low-risk groups might be partially due to different tumor immune microenvironments. To further confirm this conjecture, we conducted TMB and immune infiltration analyses.

TMB is a quintessential predictive marker for cancer immunotherapy [54]. In this study, we found that the high-risk group had higher TMB scores. Namely, patients with higher PRS might benefit more from immunotherapy than those with low PRS. We also uncovered the association between PRS and immune cell infiltration using various approaches. According to our research, the PRS was positively related to M0 and M2 macrophages and inversely associated with most other immune cells. Importantly, we utilized the “Scissor” algorithm to better illustrate the differences in immune infiltration between high and low PRS groups at the single-cell level. Concordant with former results, we found that effector memory CD8+ T cell, central memory CD4+ T cell, M1 macrophage, monocyte, B cell and plasma cell were enriched in the Scissor- (low-risk) group, while T follicular helper cell, M2 macrophage, regulatory T cell and plasmacytoid dendritic cell were more abundant in the Scissor+ (high-risk) group. Furthermore, patients in the low-risk group were more sensitive to commonly used chemotherapeutic drugs including docetaxel, epirubicin, paclitaxel, vinblastine, vincristine and vinorelbine. All of these results are in alignment with the better survival outcomes of patients in the low-risk group.

Notably, we identified CHMP2B as the hub gene among ten PANoptosis-related genes combining LASSO, WGCNA, RF and XGBoost methods. The role of CHMP2B in cancer is still controversial. On the one hand, it has a detrimental effect on the survival of patients with pancreatic adenocarcinoma [55], and overexpression of cPLA2G4A and CHMP2B simultaneously is correlated to the higher grade of myxofibrosarcoma and myxoid liposarcoma [56]. On the other hand, CHMP2B is decreased in patients with endometrial carcinoma and the urinary exosomes of patients with colorectal cancer compared to healthy individuals [57, 58]. Recently, Guo et al. revealed that cytoplasmic Yes1-associated transcriptional regulator (YAP1) could inhibit the proliferation of breast tumors by promoting autophagy, potentially through the combination of CHMP2B and VPS4B [59]. However, the experiments conducted and the model constructed in this research were mainly based on YAP1, while the role of CHMP2B itself remained to be further explored. In our study, we applied IHC analysis based on tissue microarray to further validate the role of CHMP2B in breast cancer. The results indicated that patients with highly expressed CHMP2B had significantly lower expressions of CD4 and CD8 and worse prognoses, which were consistent with bioinformatic analyses.

Despite the merits of our findings, several limitations remain to be noted. First, our model was constructed and validated retrospectively. Hence, prospective research in the real world is highly required to underpin the clinical utility of our model. Additionally, although IHC analysis based on tissue array was performed, rational and precise mechanical experiments are considered necessary to unravel the underlying mechanisms of these PANoptosis-related genes.

In conclusion, our study constructed and validated a PANoptosis-based prognostic model, which provided significant value in predicting the survival outcomes of BC patients. Besides, we processed a series of comprehensive analyses between patients with high and low PRS, further confirming the importance of PANoptosis-related gene signature in the modulation of TME and drug sensitivity in BC, providing pivotal insights for subsequent mechanical research and helping clinicians make more personalized treatment decisions.

Abbreviations

BC: breast cancer; RF: random forest; LASSO: Least absolute shrinkage and selection operator; PRS: PANoptosis-related score; CHMP2B: charged multivesicular body protein 2B; scRNA-seq: single-cell RNA sequencing; GEO: Gene Expression Omnibus; TNBC: triple-negative breast cancer; HER2+: human epidermal growth factor receptor 2-positive; ER+: estrogen receptor-positive; OS: overall survival; BN: Bayesian network; DEGs: differentially expressed genes; TMB: tumor mutation burden; IPS: Immunophenoscore; WGCNA: weighted gene co-expression network analysis; IHC: immunohistochemistry; TMA: tissue microarray; CHMP6: charged multivesicular body protein 6; TP63: tumor protein p63; IRF2: Interferon regulatory factor 2; CASP7: caspase 7; CHMP4A: charged multivesicular body protein 4A; IL-18: interleukin 18; AIFM1: apoptosis inducing factor mitochondria associated 1; GZMB: granzyme B; CHMP3: charged multivesicular body protein 3; TME: tumor microenvironment; NGS: next-generation sequencing; ACD: accidental cell death; RCD: regulated cell death; PCD: programmed cell death; KM: Kaplan-Meier; ROC: receiver operating characteristic; AUC: area under curve; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology; TCGA: The Cancer Genome Atlas; TIIC: tumor-infiltrating immune cell; ssGSEA: single sample Gene Set Enrichment Analysis; TCIA: the Cancer Immunome Atlas; TIDE: Tumor Immune Dysfunction and Exclusion; HLA: human leukocyte antigen; IC50: half-maximal inhibitory concentration; Scissor: Single-Cell Map of Diverse Immune Phenotypes in the Breast Tumor Microenvironment; GDSC: Genomics of Drug Sensitivity in Cancer; GPCR: G protein-coupled receptor; TCR: T cell antigen receptor; TOM: Topological Overlap Matrix.

Supplementary Material

Supplementary figures and tables.

Acknowledgements

Funding

This research was supported by grants from the National Natural Science Foundation of China (81872135, 82002791) and the Funds for Distinguished Young Scientists of the Second Affiliated Hospital of Harbin Medical University.

Consent to participate

Informed consent was obtained from all subjects involved in the study.

Data availability

Publicly available datasets were analyzed in this study. This data can be found here: https://xenabrowser.net/ and https://www.cbioportal.org/.

Ethics approval

This study was performed in line with the principles of the Declaration of Helsinki and approved by the Institutional Review Board of The Second Affiliated Hospital of Harbin Medical University.

Author contributions

TSY, FM and BLG designed the study. WLC, JRZ, YSL, YQD and HYZ collected the data. TSY, TW, ABH, JYF, MCL, YLL, ZYR and SSS conducted the statistical analyses and visualized the results. TSY and YHS collectively conceptualized the manuscript. FM and BLG edited the manuscript and provided critical comments. All authors read and approved the final manuscript.

Competing Interests

The authors declare that there is no conflict of interest regarding the publication of this article.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A. et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA: a cancer journal for clinicians. 2021

2. Waks AG, Winer EP. Breast Cancer Treatment: A Review. JAMA. 2019;321:288-300

3. Yu F, Quan F, Xu J, Zhang Y, Xie Y, Zhang J. et al. Breast cancer prognosis signature: linking risk stratification to disease subtypes. Brief Bioinform. 2019;20:2130-40

4. Abd-Elnaby M, Alfonse M, Roushdy M. Classification of breast cancer using microarray gene expression data: A survey. J Biomed Inform. 2021;117:103764

5. Hong M, Tao S, Zhang L, Diao L-T, Huang X, Huang S. et al. RNA sequencing: new technologies and applications in cancer research. J Hematol Oncol. 2020;13:166

6. Poorvu PD, Gelber SI, Rosenberg SM, Ruddy KJ, Tamimi RM, Collins LC. et al. Prognostic Impact of the 21-Gene Recurrence Score Assay Among Young Women With Node-Negative and Node-Positive ER-Positive/HER2-Negative Breast Cancer. J Clin Oncol. 2020;38:725-33

7. Hanahan D. Hallmarks of Cancer: New Dimensions. Cancer Discov. 2022;12:31-46

8. Galluzzi L, Vitale I, Aaronson SA, Abrams JM, Adam D, Agostinis P. et al. Molecular mechanisms of cell death: recommendations of the Nomenclature Committee on Cell Death 2018. Cell Death and Differentiation. 2018;25:486-541

9. Drago JZ, Chandarlapaty S, Jhaveri K. Targeting Apoptosis: A New Paradigm for the Treatment of Estrogen Receptor-Positive Breast Cancer. Cancer Discov. 2019;9:323-5

10. Wei X, Xie F, Zhou X, Wu Y, Yan H, Liu T. et al. Role of pyroptosis in inflammation and cancer. Cell Mol Immunol. 2022;19:971-92

11. Gong Y, Fan Z, Luo G, Yang C, Huang Q, Fan K. et al. The role of necroptosis in cancer biology and therapy. Mol Cancer. 2019;18:100

12. Malireddi RKS, Kesavardhana S, Kanneganti T-D. ZBP1 and TAK1: Master Regulators of NLRP3 Inflammasome/Pyroptosis, Apoptosis, and Necroptosis (PAN-optosis). Frontiers In Cellular and Infection Microbiology. 2019;9:406

13. Jiang W, Deng Z, Dai X, Zhao W. PANoptosis: A New Insight Into Oral Infectious Diseases. Frontiers In Immunology. 2021;12:789610

14. Christgen S, Tweedell RE, Kanneganti T-D. Programming inflammatory cell death for therapy. Pharmacol Ther. 2022;232:108010

15. Karki R, Sharma BR, Lee E, Banoth B, Malireddi RKS, Samir P. et al. Interferon regulatory factor 1 regulates PANoptosis to prevent colorectal cancer. JCI Insight. 2020 5

16. Pan H, Pan J, Li P, Gao J. Characterization of PANoptosis patterns predicts survival and immunotherapy response in gastric cancer. Clin Immunol. 2022;238:109019

17. He P, Ma Y, Wu Y, Zhou Q, Du H. Exploring PANoptosis in breast cancer based on scRNA-seq and bulk-seq. Front Endocrinol (Lausanne). 2023;14:1164930

18. Goldman MJ, Craft B, Hastie M, Repečka K, McDade F, Kamath A. et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nature Biotechnology. 2020;38:675-8

19. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO. et al. Integrative analysis of complex cancer genomics and clinical profiles using the cBioPortal. Sci Signal. 2013;6:pl1

20. Cerami E, Gao J, Dogrusoz U, Gross BE, Sumer SO, Aksoy BA. et al. The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data. Cancer Discov. 2012;2:401-4

21. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 2015;1:417-25

22. Pal B, Chen Y, Vaillant F, Capaldo BD, Joyce R, Song X. et al. A single-cell RNA expression atlas of normal, preneoplastic and tumorigenic states in the human breast. EMBO J. 2021;40:e107333

23. Liaw A, Wiener M. Classification and Regression by RandomForest. Forest. 2001 23

24. Friedman J, Hastie T, Tibshirani R. Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw. 2010 33

25. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550

26. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2:100141

27. Yu G, He Q-Y. ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Mol Biosyst. 2016;12:477-9

28. Sato N, Tamada Y, Yu G, Okuno Y. CBNplot: Bayesian network plots for enrichment analysis. Bioinformatics. 2022

29. Mayakonda A, Lin D-C, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 2018;28:1747-56

30. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7

31. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D. et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 2017;18:248-62

32. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X. et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nature Medicine. 2018;24:1550-8

33. Fu J, Li K, Zhang W, Wan C, Zhang J, Jiang P. et al. Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 2020;12:21

34. Hao Y, Hao S, Andersen-Nissen E, Mauck WM, Zheng S, Butler A. et al. Integrated analysis of multimodal single-cell data. Cell. 2021 184

35. Hu C, Li T, Xu Y, Zhang X, Li F, Bai J. et al. CellMarker 2.0: an updated database of manually curated cell markers in human/mouse and web tools based on scRNA-seq data. Nucleic Acids Research. 2023;51:D870-D6

36. Azizi E, Carr AJ, Plitas G, Cornish AE, Konopacki C, Prabhakaran S. et al. Single-Cell Map of Diverse Immune Phenotypes in the Breast Tumor Microenvironment. Cell. 2018 174

37. Sun D, Guan X, Moran AE, Wu L-Y, Qian DZ, Schedin P. et al. Identifying phenotype-associated subpopulations by integrating bulk and single-cell sequencing data. Nature Biotechnology. 2022;40:527-38

38. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S. et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Research. 2013;41:D955-D61

39. Wishart DS, Feunang YD, Guo AC, Lo EJ, Marcu A, Grant JR. et al. DrugBank 5.0: a major update to the DrugBank database for 2018. Nucleic Acids Research. 2018;46:D1074-D82

40. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559

41. Chen T, Guestrin C. Xgboost: A scalable tree boosting system. Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining. 2016 p. 785-94

42. Camp RL, Dolled-Filhart M, Rimm DL. X-tile: a new bio-informatics tool for biomarker assessment and outcome-based cut-point optimization. Clinical Cancer Research: an Official Journal of the American Association For Cancer Research. 2004;10:7252-9

43. Tian W, Luo Y, Tang Y, Kong Y, Wu L, Zheng S. et al. Novel Implication of the Basement Membrane for Breast Cancer Outcome and Immune Infiltration. Int J Biol Sci. 2023;19:1645-63

44. Cui H, Ren X, Dai L, Chang L, Liu D, Zhai Z. et al. Comprehensive analysis of nicotinamide metabolism-related signature for predicting prognosis and immunotherapy response in breast cancer. Frontiers In Immunology. 2023;14:1145552

45. Li W, Wang X, Li C, Chen T, Zhou X, Li Z. et al. Identification and validation of an m6A-related gene signature to predict prognosis and evaluate immune features of breast cancer. Hum Cell. 2023;36:393-408

46. Bedoui S, Herold MJ, Strasser A. Emerging connectivity of programmed cell death pathways and its physiological implications. Nat Rev Mol Cell Biol. 2020;21:678-95

47. Strasser A, Harris AW, Huang DC, Krammer PH, Cory S. Bcl-2 and Fas/APO-1 regulate distinct pathways to lymphocyte apoptosis. EMBO J. 1995;14:6136-47

48. Cho YS, Challa S, Moquin D, Genga R, Ray TD, Guildford M. et al. Phosphorylation-driven assembly of the RIP1-RIP3 complex regulates programmed necrosis and virus-induced inflammation. Cell. 2009;137:1112-23

49. Cai Z, Jitkaew S, Zhao J, Chiang H-C, Choksi S, Liu J. et al. Plasma membrane translocation of trimerized MLKL protein is required for TNF-induced necroptosis. Nat Cell Biol. 2014;16:55-65

50. Christgen S, Kanneganti T-D. Inflammasomes and the fine line between defense and disease. Curr Opin Immunol. 2020;62:39-44

51. Maitra S, Bhattacharya D, Paul S, Ghosh Chowdhury P, Mandal D, Haldar PK. et al. Programmed Cell Death Protein 1 (PD-1) in relation to PANoptosis: Immune Pharmacological Targets for Management of Breast Adenocarcinoma. Endocr Metab Immune Disord Drug Targets. 2023

52. Gatti V, Bongiorno-Borbone L, Fierro C, Annicchiarico-Petruzzelli M, Melino G, Peschiaroli A. p63 at the Crossroads between Stemness and Metastasis in Breast Cancer. International Journal of Molecular Sciences. 2019 20

53. Zhou T, Damsky W, Weizman O-E, McGeary MK, Hartmann KP, Rosen CE. et al. IL-18BP is a secreted immune checkpoint and barrier to IL-18 immunotherapy. Nature. 2020;583:609-14

54. Chan TA, Yarchoan M, Jaffee E, Swanton C, Quezada SA, Stenzinger A. et al. Development of tumor mutation burden as an immunotherapy biomarker: utility for the oncology clinic. Ann Oncol. 2019;30:44-56

55. Zhang Q, Lv L, Ma P, Zhang Y, Deng J, Zhang Y. Identification of an Autophagy-Related Pair Signature for Predicting Prognoses and Immune Activity in Pancreatic Adenocarcinoma. Frontiers In Immunology. 2021;12:743938

56. Wang T, Goodman MA, McGough RL, Weiss KR, Rao UNM. Immunohistochemical analysis of expressions of RB1, CDK4, HSP90, cPLA2G4A, and CHMP2B is helpful in distinction between myxofibrosarcoma and myxoid liposarcoma. Int J Surg Pathol. 2014;22:589-99

57. Ma L, Yu H, Zhu Y, Xu K, Zhao A, Ding L. et al. Isolation and proteomic profiling of urinary exosomes from patients with colorectal cancer. Proteome Sci. 2023;21:3

58. Yang Y, Wang M. Genomic analysis of the endosomal sorting required for transport complex III pathway genes as therapeutic and prognostic biomarkers for endometrial carcinoma. Transl Cancer Res. 2022;11:3108-27

59. Guo Y, Cui Y, Li Y, Jin X, Wang D, Lei M. et al. Cytoplasmic YAP1-mediated ESCRT-III assembly promotes autophagic cell death and is ubiquitinated by NEDD4L in breast cancer. Cancer Commun (Lond). 2023;43:582-612

Author contact

Corresponding authors: Baoliang Guo (baoliangguo2020edu.cn); Fei Ma (wafsfdcom).

Received 2023-9-13
Accepted 2023-11-14
Published 2024-1-1

Citation styles

©2024 Ivyspring International Publisher. Terms of use