Consensus clustering analysis of PANoptosis-related genes

A consensus clustering analysis based on PRGs expression was performed to investigate the connections between PRGs and colon cancer subtypes using the ConsensusClusterPlus R package. The classification with the highest intragroup correlations and the lowest intergroup correlations was identified. The prognosis of the two clusters was compared using the Kaplan-Meier (KM) method and log-rank test. Principal component analysis (PCA) was performed using the stats R package. Differences in clinical features between two clusters were analyzed using the Wilcoxon test, and DEGs between two clusters were identified with the criteria of |log fold change (FC)| >1 and p-value < 0.05, using limma package. To explore the differences in biological processes between the two PRG clusters, we performed gene set variation analysis (GSVA) using gsva R package. Single-sample gene set enrichment analysis (ssGSEA) was used to calculate the scores of infiltrating immune cells and evaluate the activity of immune-related functions.

Gene oncology and kyoto encyclopedia of genes and genomes analyses

To understand the biological functions and pathways related to the DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the ggplot2, Bioconductor, and org. Hs.eg.db R packages. p-values and q-values < 0.05 were considered statistically significant.

Construction of the PANoptosis-related prognostic signature

Prognosis-related DEGs (PRDEGs) were selected using a univariate Cox regression analysis. To identify additional PANoptosis-related genes for signature construction. We classified patients into two distinct clusters based on the expression of PRDEGs. Survival time, clinical features, and PRG expression were compared between the two gene clusters, and DEGs between gene clusters A and B were identified. Finally, seven genes were included to construct the prognostic signature after least absolute shrinkage and selection operator (LASSO) regression analysis and multivariate Cox regression analysis using the survival, survminer, and glmnet R packages. The risk score was calculated based on the expression levels of the seven genes, and the patients were divided into high- and low-risk groups using the median risk score. The KM analysis was used to evaluate survival differences between high-risk and low-risk groups, and the receiver operating characteristic (ROC) and area under the curve (AUC) were used to test the prediction efficiency of the risk score. A nomogram model was developed based on risk scores and clinical factors. Calibration graphs were constructed to show the differences between the nomogram-predicted and actual survival rates of colon cancer patients.

Evaluation of the tumor microenvironment between the high- and low-risk groups

To better understand the relationship between the risk score and tumor microenvironment (TME), CIBERSORT was used to quantify the abundance of infiltrating immune cells in heterogeneous samples from the high- and low-risk groups. The correlation between the abundance of immune cells and the seven genes was analyzed. TME scores, including stromal, immune, and ESTIMATE scores, of high- and low-risk groups were also compared using the Wilcoxon signed-rank test.

Analyses of mutations, microsatellite instability and cancer stem cell index between high- and low-risk groups

To explore the gene mutations in colon cancer patients in high- and low-risk groups, we generated the mutation annotation format (MAF) using the maftools R package. The tumor mutation burden (TMB) score of patients was calculated, and the correlation between the risk score and TMB was analyzed using the Spearman method. Furthermore, we analyzed the association between risk groups and MSI and CSC index.

Response to immunotherapy and chemotherapeutic drugs

The Cancer Immunome Atlas (TCIA, https://tcia.at/) is a dataset containing TCGA data of 20 solid cancers and more than 8,000 tumor samples. It can be used to calculate the immunophenotypic score (IPS) of tumor samples to predict the response to cytotoxic T lymphocyte antigen-4 (CTLA-4) and programmed cell death protein 1 (PD-1) blockers (Charoentong et al., 2017). IC50 is half of the maximum inhibitory concentration and represents the concentration of the drug required for 50% inhibition of cancer cells. To determine the relationship between the risk score and response to immunotherapy and chemotherapeutic drugs, checkpoint expression, immune subtypes, IPS of tumor samples, and the IC50 of drugs in the two risk groups were calculated and compared.

Identification of PANoptosis-related gene clusters in colon cancer

To explore the interactions between the 19 PRGs and their prognostic significance, a network was constructed, as shown in Figure 2A. Kaplan–Meier curves of the relationship between PRGs expression and the prognosis of colon cancer patients were shown in Supplementary Figure S1. Consensus clustering analysis was performed to explore the relationship between PRG expression and tumor classification (Supplementary Figure S2). Clusters with the highest intragroup correlations and lowest intergroup correlations were identified. By increasing the clustering variable (k), we found that when k = 2, classification met the standard. Colon cancer patients were divided into two PRG clusters (A and B) based on PRG expression levels (Figure 2B). As shown in Figure 2C, patients in PRG cluster A had a significantly longer survival time than those in cluster B (p = 0.048). PCA showed a satisfactory separation between PRG cluster A and B (Figure 2D). Figure 2E shows the association between PRG clusters and clinical features and PRG expression in colon cancer patients. Tumor infiltration and lymph node metastasis correlated with PRG clusters (p < 0.05). GSVA showed that PRG cluster A was significantly enriched in immune-related pathways, including natural killer cell-mediated cytotoxicity, antigen processing and presentation, primary immunodeficiency, B cell, and T-cell receptor signaling pathways (Figure 2F). To evaluate the differences in immune cell infiltration between the two clusters, ssGSEA was performed, and the results showed that PRG cluster A had higher immune cell infiltration levels, including those of activated B cells, activated CD4 + T-cells, activated CD8 + T-cells, activated dendritic cells, macrophages, mast cells, and natural killer cells (Figure 2G).

Open in a separate window

FIGURE 2

PRGclusters and clinical characteristics between colon cancer samples in two clusters. Relationship of tumor microenvironment in two PRGclusters. (A) Interactions among PRGs in colon cancer. The lines among the PRGs represents their interactions. Blue and red represent negative and positive correlations. (B) Two PRGclusters were defined using consensus clustering analyses. (C) KM curve indicated that PRGcluster a had longer survival time than PRGcluster B (p = 0.048). (D) PCA showed good distiction between two PRGclusters. (E) Heatmaps showed the relationship between PRGclusters and clinical features and PRGs expression in colon cancer patients. (F) GSVA showed the enriched pathways in PRGclusters. (G) ssGSEA investigated the differences of immune cell infiltration between two clusters. *p < 0.05; **p < 0.01.

Identification of gene clusters based on differentially expressed genes

DEGs were identified, and GO and KEGG analyses showed the relevant biological processes (BP), cellular components (CC), molecular functions (MF), and pathways (Figures 3A,B). These DEGs were mainly related to the BP of T cell activation, leukocyte cell-cell adhesion, and response to interferon-gamma, and were correlated with the CC such as the external side of the plasma membrane, major histocompatibility complex (MHC) class II protein complex, and MHC protein complex. Furthermore, they were involved in the MF of immune receptor activity, chemokine activity, and antigen binding. According to KEGG analysis, these DEGs participate in certain cancer-related pathways, including the chemokine signaling pathway, NOD-like receptor signaling pathway, and NF-κB signaling pathway. PRDEGs were identified using univariate Cox regression analysis. Patients were then divided into two clusters (gene cluster A and gene cluster B) based on PRDEG expression (Supplementary Figure S3; Figure 3C) shows that cluster A had higher survival rates than cluster B (p = 0.002). In Figure 3D, the boxplot shows that FADD, CASP6, CASP7, IRF1, AIM2, ZBP1, CASP1, RIPK1, RIPK3, and TRADD were upregulated in cluster A, whereas TAB3 and PARP1 were downregulated in cluster B (p < 0.05). The heatmaps show the association between gene clusters and clinical features and between PRDEGs expression and PRG clusters. Furthermore, gene clusters were significantly related to tumor infiltration and metastasis (p < 0.05) (Figure 3E).

Open in a separate window

FIGURE 3

Identification of geneclusters based on DEGs. (A–B) GO and KEGG analyses showed the relevant biological processes (BP), cellular components (CC), molecular functions (MF) and pathways. (C) Heatmap showed the association between genecluster and clinical features. (D) KM curves showed genecluster A had a more favorable prognosis. (E) Expression levels of PRGs in two geneclusters. *p < 0.05; **p < 0.01; and ***p < 0.001.

Development abd validation of the PANoptosis-related prognostic signature

LASSO and Cox regression analyses were performed to screen for prognosis-related DEPRGs (Figures 4A,B). After selection, seven genes were included in the calculation of the risk score based on the following formula: Risk score = ${\displaystyle \sum }_{i=1}^{\mathrm{n}}\beta i*\lambda i$ , where n represents the number of genes included to construct the signature, and $\beta i$ and $\lambda i$ represent the regression coefficient and gene expression value, respectively. Boxplots showed that PRG cluster B and gene cluster B had higher risk scores than PRG cluster A and gene cluster A (Figures 4C,D). A Sankey diagram showed the associations among PRG cluster, gene cluster, risk groups, and survival status (Figure 4E). Fifteen of the nineteen PRGs were differentially expressed between the high- and low-risk groups (Figure 4F). Seven genes were used to construct the prognostic signature; Figure 5A shows the expression differences of these seven genes between the two risk groups. Based on the risk score, patients were divided into high- and low-risk groups, and patients with higher risk scores had a higher risk of mortality (Figure 5B). As shown in Figure 5C, the KM curve was plotted to show the survival differences between the two groups. Patients in the high-risk group had a significantly lower probability of survival than those in the low-risk group (p < 0.001). ROC curves were drawn to test the prediction efficiency of the risk score, and the AUCs for one-, three-, and 5-years survival were 0.612, 0.650, and 0.676, respectively (Figure 5D). The results of the risk score in TCGA (Supplementary Figure S4) and GSE39582 (Supplementary Figure S5) cohorts were also shown. The risk score and other clinical features were used to construct a nomogram model (Figure 5E). Calibration plots showing the differences between the nomogram-predicted and actual survival probabilities of colon cancer patients showed that the predicted survival probabilities were close to the actual survival probabilities (Figure 5F), indicating that this nomogram model accurately predicted the survival of colon cancer patients. The results for three independent validation cohorts, which were GSE17536 (Figure 5G, p = 0.041, 1-year AUC = 0.598, 3-years AUC = 0.624, 5-years AUC = 0.589), GSE17537 (Figure 5H, p = 0.048, 1-year AUC = 0.728, 3-years AUC = 0.624, 5-years AUC = 0.542), and GSE29621 (Figure 5I, p = 0.011, 1-year AUC = 0.763, 3-years AUC = 0.717, 5-years AUC = 0.702), revealed that the risk score could efficiently predict patient survival.

Open in a separate window

FIGURE 4

Identification of 7 genes for calculating the risk score and the relationship between molecular classifications, PRG expression levels and the risk score. (A–B) The LASSO regression analysis and partial likelihood deviance on the prognostic genes. (C–D) Association between risk score and molecular classifications. (E) Sankey plot showed the correlation between molecular classifications, risk groups and survival status in colon cancer patients. (F) Expression levels of PRGs in two risk groups. *p < 0.05; **p < 0.01; and ***p < 0.001.

Open in a separate window

FIGURE 5

Construction and validation of the prognostic signature. (A) Heatmap showed the expression of 7 genes in two risk groups. (B) Risk score and survival outcome of each case. (C) KM curve showed that patients in high-risk group had a worse prognosis. (D) The AUC for 1-, 3- and 5-years survival were 0.612, 0.650, and 0.676, respectively. (E) Nomogram using risk score and other clinical features were constructed for predicting survival of colon cancer patients. (F) Calibration graphs investigated that the actual survival rates of colon cancer patients were close to the nomogram-predicted survival rates. The KM and ROC methods were used to evaluate the efficiency of the risk score at predicting patient survival in GSE17536 (G), GSE29621 (H), and GSE38832 (I) CRC datasets.

Comparative evaluation of the tumor microenvironments of high- and low-risk groups

Figure 6Ashows the correlation between the risk score and immune cell abundance: M0 macrophages, M2 macrophages, activated mast cells, and neutrophils were positively related to the risk score, while naive B cells, activated dendritic cells, resting dendritic cells, M1 macrophages, resting mast cells, resting NK cells, plasma cells, activated memory CD4 + T-cells, resting memory CD4 + T-cells, CD8 + T-cells, and follicular helper T-cells were negatively related to the risk score. The relationship between the abundance of immune cells and seven genes in the prognostic signature was also evaluated (Figure 6B). The TME scores in the two risk groups were calculated, and the high-risk group was found to have a higher stromal and lower immune score (Figure 6C).

Open in a separate window

FIGURE 6

Evaluation of tumor microenvironment in high- and low-risk groups. (A) Relationship between risk score and different immune cell types. (B) Correlation between the abundance of immune cells and seven genes in the prognostic signature. (C) Correlation between risk score and immune-related scores. *p < 0.05; **p < 0.01; and ***p < 0.001.

Comparative analysis of mutations, microsatellite instability and cancer stem cell index in high- and low-risk groups

Differences in somatic mutations between the two risk groups of colon cancer patients were analyzed; the five most mutated genes in the high- and low-risk groups were APC, TP53, TTN, KRAS, and SYNE1 (Figures 7A,B). TMB (Figure 7C) and MSI (Figure 7D) did not show a significant relationship with the risk score, while CSC (Figure 7E) was negatively correlated with the risk score (R = −0.15, p < 0.01).

Open in a separate window

FIGURE 7

Comprehensive analyses of risk score in colon cancer. The somatic gene mutations in high-risk group (A) and low-risk group (B). TMB (C) and MSI (D) did not show significant correlations with risk score while CSC (E) negatively correlated with risk score.

Response to immunotherapy and chemotherapeutic drugs

To analyze the ability of risk score to predict potential checkpoint blockade therapy, boxplots were drawn to show the differences in immune checkpoint gene expression between the high- and low-risk groups (Figure 8A). Checkpoint genes, including CTLA4, LAG3, ID O 2, CD274, and PDCD1, had higher expression levels in low-risk groups. In Figure 8B, cluster1 (C1), C2, C3, and C4 represent wound healing, IFN-gamma dominant, inflammatory, and lymphocyte depleted immune subgroups, respectively (Thorsson et al., 2019). The results showed that C3 samples were almost equally distributed between the two groups, but there were more C1 and C4 samples and fewer C2 samples in the high-risk subgroup than in the low-risk subgroup. Violin plots showed the relationship between IPSs and risk groups; a higher IPS represented a better response to PD-1 and CTLA-4 blockers (Figure 8C). We also found that eight drugs had lower IC50 values in the high-risk group, including bexarotene, bicalutamide, dasatinib, doxetacel, elesclomol, imatinib, midostaurin, and pazopanib (Figure 8D).

Open in a separate window

FIGURE 8

Response to anti-tumor tharapy of colon cancer patients in two risk groups. (A) The differences of immune checkpoint gene expression in high-risk and low-risk groups. (B) Heatmap and table showing the distribution of colon cancer immune subtypes between two risk groups. (C) Violin plots showed the relationship between IPSs and risk groups. (D) Eight therapeutic drugs showed significant IC50 differences. *p < 0.05; **p < 0.01; and ***p < 0.001.

We acknowledged TCGA and GEO database for providing their platform and contributors for uploading their meaningful datasets.