We initiated our research by acquiring a comprehensive dataset of transcriptome sequencing and clinical data from The Cancer Genome Atlas-Colorectal Cancer (TCGA-COAD) repository. The dataset included 458 unique patient cases, comprising 41 normal tissue samples and 483 COAD tissue samples. Strict quality control was implemented to ensure the reliability and consistency of the procured samples, excluding the datasets with incomplete records or ambiguous clinical outcomes. For a thorough investigation of PANoptosis, we incorporated five distinct pathways from the Molecular Signatures Database (MSigDB) version 7.4 (GSEA | MSigDB: gsea-msigdb.org). These pathways included REACTOME_PYROPTOSIS, HALLMARK_APOPTOSIS, KEGG_APOPTOSIS REACTOME_APOPTOSIS, and KEGG_NECROPTOSIS ( Supplementary Table S1 ). Additionally, based on earlier studies, we identified a set of immune-related genes that are pivotal to our analysis (28).

The interaction between PANoptosis and the immune system in cancer has been previously reported (29, 30). To obtain a clearer understanding of this interaction, the transcriptome data were analyzed using single-sample gene set enrichment analysis (ssGSEA) to generate enrichment scores for PANoptosis-related and immune-related genes. This process was conducted using the GSVA package in R (31). ssGSEA is a widely used method to assess gene set enrichment within individual samples, which allows the identification of specific molecular signatures.

To investigate the potential association between PANoptosis and immune response, we performed a correlation analysis between enrichment scores of PANoptosis-related and immune-related genes using the ‘cor.test’ function in R. A strong correlation between immune features and PANoptosis might suggest that these processes are interconnected in initiation and progression of COAD. For the strongly correlated features (Spearman correlation coefficient greater than 0.7), we utilized hierarchical clustering using the ‘hclust’ function in R (32). This technique groups samples exhibiting similar patterns, enabling the assessment of distinct molecular subtypes. Strong correlations reveal co-expression between immune response and PANoptosis-related genes, while clustering identifies co-expressed gene groups.

Characterizing the differences in gene expression between identified clusters can provide vital insights into the underlying biological pathways. Differential gene expression analysis between the clusters was performed using the ‘limma’ package in R (33), with a cutoff of logFC = 0.5 and p-value < 0.05.

The survival analysis for PANoptosis genes was performed using the ‘survival’ and ‘survminer’ packages in R. By identifying genes associated with survival differences, we can gain valuable information on potential prognostic biomarkers. The Cox proportional hazards model was applied for univariate survival analyses. This analysis enable the investigation of the relationship between gene expression levels and survival outcomes, providing essential information for prognostic assessment.

We aimed to identify the genes that are both differentially expressed in the immune-related groups and associated with survival. The intersection of differential genes and PANoptosis-related prognostic genes was calculated. To further explore the molecular subtypes for better clinical assessment, we used consensus clustering using the ‘ConsensusClusterPlus’ package in R with the intersected gene set. The maximum number of clusters (K) to test was set at 9. A total of 50 resampling iterations were performed for each K. In each iteration, 80% of the samples (pItem=0.8) were randomly selected.

Tumor mutation burden (TMB) is a measure of the number of somatic mutations present in a tumor. Higher TMB levels correlate with a better immunotherapy response, making it a valuable metric for assessing treatment potential. The TMB for each sample group was calculated using the maftools package in R based on the mutation data. The ‘ESTIMATE’ package in R was used to infer the stromal and immune scores, reflecting the immune landscape. Understanding the immune landscape within tumors allows for improved tumor characterization and prediction of response to immunotherapies.

Tumor immune escape is a significant concern when considering immunotherapy treatment strategies. To predict the response to immune checkpoint blockades, we utilized the web tool Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/), which calculates immune evasion scores based on gene expression profiles.

To explore the cellular context of genes constituting the identified biomarker signature, single-cell transcriptomic data was analyzed using The Tumor Immune Single-Cell Hub (TISCH, http://tisch.comp-genomics.org). By investigating the gene signature expression in different cell types, we can extract valuable biological insights into how these genes may functionally affect the tumor microenvironment.

Constructing a prognostic model for clinical applications is critical for patient stratification and treatment decision-making. The Least Absolute Shrinkage and Selection Operator (LASSO) method was utilized to minimize overfitting during the establishment of the prognostic model using the ‘glmnet’ package in R (34). Additionally, we employed Cox regression to evaluate relationships between survival time and predictors. This approach provides a robust and interpretable model, allowing clinicians to evaluate patients’ risk profiles.

We assessed the predictive accuracy of our risk score model using the ROC curve (‘pROC’ package in R) and Calibration curve (‘rms’ package in R) analyses (35). These methods offer quantitative measurement of our proposed model’s performance and allow comparison with other clinical factors.

To enhance our understanding of the biological processes and pathways altered within the identified gene sets, differentially expressed genes were subjected to both GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analyses using the clusterProfiler package in R (36). These approaches provide a systematic representation of the biological functions modulated by the genes of interest.

A nomogram was constructed, integrating the risk score and other relevant clinicopathological factors, to provide a comprehensive representation for the prediction of individualized survival probability. The nomogram was built using the ‘rms’ package in R. Moreover, the performance of the nomogram was further assessed by plotting calibration curves. Calibration curves offer a visual comparison between the predicted and observed outcomes, providing a direct evaluation of the accuracy of the nomogram.

The study flowchart is depicted in Figure 1 . In this study, we first performed a comprehensive Single-Sample Gene Set Enrichment Analysis (ssGSEA) for the PANoptosis gene set and immune-related gene set, providing an integrated view of the PANoptosis and immune landscape in the COAD sample set. We observed a strong correlation (Spearman correlation coefficient > 0.7) between the enrichment scores of the PANoptosis gene set and four immune feature gene sets (Treg, parainflammation, CCR, and immune checkpoint), as depicted in Figure 2A . The strong correlation indicates a likely co-regulation or mutual influence. Immune responses are intricately tied to PANoptosis, a process that maintains cellular homeostasis. For instance, Tregs, known to suppress immune responses, might be leveraging PANoptosis to maintain immune tolerance (37). Parainflammation, a response to tissue stress, could be interlinked with PANoptosis as a means of removing damaged cells and promoting tissue repair (38). CCR and immune checkpoint genes, key mediators of immune responses, might be influencing or influenced by apoptotic processes, contributing to the immune landscape (39, 40). However, the precise mechanisms remain to be elucidated and warrant further investigation.

Upon identifying immune feature gene sets with a Spearman correlation coefficient greater than 0.7, hierarchical clustering was performed. This resulted in the categorization of 189 samples into cluster 1 and 294 samples into cluster 2. The distribution of the enrichment scores for the four immune-related gene sets - Treg, Parainflammation, CCR, and checkpoint - across the two clusters is depicted in Figure 2B . Notably, all these scores were significantly higher in cluster 2, implying a heightened immune response or activity in this cluster. The clear separation into two clusters suggests inherent molecular differences between them.

A differential gene expression analysis was performed between the two clusters with selection criteria of logFC= 0.5 and p-value< 0.05. As visualized in the volcano plot ( Figure 2C ), 2307 genes exhibited higher expression in cluster 1 (characterized by high immune feature gene set enrichment scores), while 741 genes were upregulated in cluster 2. The greater number of differentially highly expressed genes in cluster 1 can be attributed to the pronounced immune activity and potential interplay with various cellular pathways. A survival analysis of the PANoptosis gene set was executed using single-gene Cox analysis. The forest plot, depicted in Figure 2D , highlights the significance of certain PANoptosis genes in relation to patient survival. The intersection between the differential gene set and the survival-related PANoptosis genes yielded 10 genes of interest: PLCB2, CAV1, DAPK1, GPX3, IGFBP6, TIMP1, PMAIP1, GADD45B, ENO2, and PYGL. In light of previous scientific inquiries, it’s compelling to mention that many of the genes within this list of 10, specifically CAV1, DAPK1, GPX3, IGFBP6, TIMP1, GADD45B and ENO2 have been found associated intimately with colorectal cancer, as documented by several research studies (41–47). These genes, thereby, are candidates for further investigation concerning their role in COAD progression and therapeutic potential. This intersection is graphically represented in Figure 3A .

Based on the intersection genes identified in the previous step, consensus clustering was conducted. Figure 3B showcases the clustering results, and it was observed that the optimal classification was achieved when k=2. Consequently, 181 samples were allocated to cluster A, while 277 samples were designated to cluster B. It is worth noting that these two distinct territories, or clusters, may represent two different biological subtypes within the vast field of colorectal cancer respectively. Each subtype, potentially, could have unique implications for disease progression and prognosis. This result could enhance our molecular knowledge of colorectal cancer, potentially leading to more precise disease classification and individualized treatment strategies.

Survival analysis, a cornerstone in oncological research, was employed to fathom the prognostic implications of our identified clusters. The results, illustrated in Figures 3C (Overall Survival) and 3D (Progression-free Survival), clearly indicate that cluster B may represent a biological subtype associated with lower disease progression risk and better survival outcomes. Such a discovery could provide important insights for clinical decision-making and patient management.

To delve deeper into the genetic landscape of the tumors, we embarked on an analysis of the tumor mutation burden (TMB) across the two clusters. TMB, a measure of the number of mutations within a tumor genome, has been spotlighted for its potential prognostic value, especially in the realm of immunotherapy. Figures 4A , B offer a comparative view of the TMB across clusters A and B. Intriguingly, cluster A, as accentuated in Figure 4C , boasted a higher TMB. Additionally, we observed that cluster B had a higher frequency of APC and TP53 mutations compared to cluster A. Conversely, cluster A had a higher frequency of TTN and RYR3 mutations compared to cluster B. These findings suggest that these two subtypes seem to have different genetic mutations, which could have unique implications for disease progression and prognosis. Based on the distinct mutation frequencies of genes, we can take specific treatments for different subtypes of COAD, improve diagnostics, hone immunotherapy precision and refine clinical trial designs.

An in-depth analysis of the immune microenvironment of the two clusters was conducted. Figure 4D reveals that cluster A had elevated stromalScore, immuneScore, and ESTIMATEScore compared to cluster B, indicating a more pronounced immune and stromal activity in cluster A. The ESTIMATEScore, a combination of the stromalScore and immuneScore, gives an overall impression of the tumor purity. The reason for conducting such an analysis is to gain insight into the potential mechanisms of immune response within these clusters. The elevated values in cluster A may suggest a more reactive immune environment, potentially leading to a more effective immune response against tumor cells. The results interestingly hint at cluster B having a better prognosis. This could be surprising at first glance, as a higher level of immune and stromal activity in cluster A could be interpreted as a more robust defense against cancer progression. However, this might not always be the case. The tumor microenvironment is a complex system where a high level of immune activity can sometimes be associated with immune evasion or suppression tactics used by cancer cells. Therefore, it seems that the less pronounced immune and stromal activity in cluster B might result in a better balance between tumor-suppressive and tumor-promoting activities, leading to a more favorable prognosis. This hypothesis aligns with the notion of ‘immunoediting’, a process that shapes the immunogenicity of the tumor (48). Further investigations are needed to validate these results and explore the underlying mechanisms that contribute to the apparent paradox. Additionally, in Figure 4E , our findings, which align with the observations above, show that in cluster A, the majority of these immune function scores were indeed upregulated. To further our understanding of the tumor’s interplay with the immune system, we turned our attention to the potential for immune escape. The Tumor Immune Dysfunction and Exclusion (TIDE) framework was employed for this purpose. Consistent with the anticipated outcomes from our preceding analysis, Figure 4F , showcasing our findings, revealed a heightened TIDE score for cluster A, implying a pronounced potential for immune escape.

A single-cell transcriptomics analysis was conducted to delve deeper into the expression patterns of the 10-gene signatures. Interestingly, the 10-gene signature was observed to be highly expressed in fibroblast cells, as illustrated in Figures 5A–H . The high expression of these genes in fibroblast cells suggests that fibroblast cells may not only participate in the construction of the tumor microenvironment, but also play a role in regulating immune actions and cell death. It is intriguing, the recent studies indicate that fibroblasts hold a significant role in immune-related actions including immune escape and the process of PANoptosis, shedding light on the intricate connection of cellular mechanisms (49, 50). Additionally, the high expression of these genes may affect how these cells interact with tumor cells and other immune cells, further affecting the progression of colon cancer.

A differential analysis was undertaken between clusters A and B with a stringent selection criterion of logfc=0.5 and p-value<0.01. The volcano plot in Figure 6A reveals that 292 genes were downregulated in cluster A, while 1746 genes were upregulated.

To further our understanding of the biological underpinnings of our differential gene set, we embarked on a comprehensive enrichment analysis. This entailed both Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses ( Figures 6B, C ). The intent was to categorize our genes into meaningful biological processes, molecular functions, cellular components, and pathways. According to the KEGG enrichment analysis, the following pathways were found to be significantly enriched: Cytokine-cytokine receptor interaction, Cell adhesion molecules, Chemokine signaling pathway, MAPK signaling pathway, Focal adhesion, and Calcium signaling pathway. These three pathways, namely, Cytokine-Cytokine Receptor Interaction, Chemokine Signaling Pathway, and Cell Adhesion Molecules (CAMs), are linked with immune activities, while the MAPK signaling pathway and the Cytokine-cytokine receptor interaction are closely related to PANoptosis. Additionally, it is interesting to note that many of these pathways are related to fibroblasts, which is consistent with our previous results. The Gene Ontology (GO) enrichment analysis revealed significant pathways across three main ontology categories. For Biological Processes, pathways such as extracellular matrix organization, extracellular structure organization, and leukocyte migration were predominant. In the Cellular Component category, collagen-containing extracellular matrix, collagen trimer, and membrane raft were highlighted. Meanwhile, the Molecular Function category emphasized pathways like extracellular matrix structural constituent and glycosaminoglycan binding. These findings provide key insights into the genetic basis of our study and set the stage for deeper investigations.

Using the results from Lasso regression and Cox regression, a robust prognostic model was constructed. Figures 7A , B , representing the pathway diagram and deviation diagram respectively, aided in the optimal selection of the regularization strength λ (Lambda) value, which was determined to be 10. Finally, a risk score was constructed based on the expression levels of six genes: risk score=MAPK12 expression×0.375+ATP6V1C2 expression×0.231+HOXC11 expression×0.166+HOXD9 expression×0.224+TRPM5 expression×0.417+EEF1A2 expression×0.270.

As depicted in Figures 7C–E , a correspondence is observed between the increment in risk scores and the reduction in patients’ survival time. Six genes, namely MAPK12, ATP6V1C2, HOXC11, HOXD9, TRPM5, and EEF1A2, show elevated expression in the high-risk group. The analysis shows that the risk score is a significant predictor of patient survival time. These genes may play a critical role in the development and progression of the disease and could be potential targets for therapeutic intervention.

The ROC curve analysis for the risk score was conducted for the training set, validation set, and the entire sample set, as visualized in Figures 8A–C . In the entire group, 1-year AUC = 0.732, 3-year AUC = 0.707, and 5-year AUC = 0.729. In the present investigation, all AUC indicators surpass the benchmark of 0.7. This intimates a satisfactory discriminatory capacity of the risk score, providing a reasonably robust forecast of risk for future periods of 1 year, 3 years, and 5 years respectively. Additionally, Figure 8D presents the ROC curve correlating the risk score with clinical features, implying the risk score is a reliable predictor of future outcomes. Moreover, both univariate and multivariate Cox analyses were executed for the risk score. The forest plots in Figures 8E , F clearly demonstrate that the p-values for both analyses were less than 0.001, signifying the robustness of the risk score as a predictor.

To gauge the consistency and reliability of our risk score, we utilized the Consistency Index (C-index). Figure 8G , showcasing our findings, revealed a high C-index for our risk score, hinting at its reliability and accuracy.

In the final part of our analysis, we embarked on constructing a nomogram. This graphical tool offers clinicians a user-friendly interface to predict patient outcomes based on various factors. Figures 8H , I visually present the nomogram and its accompanying calibration curve. The calibration curve, a testament to the nomogram’s predictive accuracy with a C-index of 0.774, also further validates the efficacy of our prognostic model.